I am working on an optimization problem but I am not sure if the problem can be formulated as an integer programming problem.
Assume the cost minimization problem for a set of subscribers and providers $\mathcal{N} =\{ 1,\dots, h, \dots, i,\dots, N \}$ $$ \sum_{i=1}^{N} \sum_{j=1}^{N} m_{ij}.c_{ij} $$ in which $c_{ij} \in \mathbf{C}$ and $m_{ij} \in \mathbf{M}$.
$\mathbf{M}$ is a binary $N\times N $ matrix such that $m_{ij}=1$ says subscriber $i$ is subscribed to service provider $j$ ans its cost $c_{ij}$ is defined as
$$ c_{ij} = \frac{a_{ij}}{1+\sum_{n=1}^{\color{red}{i}} m_{nj}} - \frac{ m_{hj}. a_{hj}}{1+\sum_{n=1}^{h} m_{nj}}$$ in which $a_{ij}$ and $a_{hj}$ are given. In other words, the cost of subscribing to provider $j$ for subscriber $i$ depends on how many subscribers have already been subscribed to $j$ before $i$, and if a particular subscriber $h$ is subscribed to $j$.
In fact, such a matrix has to be found $$ \begin{bmatrix} m_{11} & \dots & m_{1j}& \dots& m_{1N} \\ \vdots & \dots & \vdots & \dots & \vdots \\ m_{i1} & \dots & m_{ij} & \dots& m_{iN} \\ \vdots & \dots & \vdots & \dots & \vdots \\ m_{N1}& \dots & m_{Nj}& \dots& m_{NN} \\ \end{bmatrix}. $$
I wondered if the optimum result can be found by integer programming or any other approaches.
Update 1: What if we have a more complicated $c_{ij}$. A $c^\prime_{ij}$ like
$$ c^\prime_{ij} = \frac{ a_{ij}}{1+\sum_{n=1}^{i} m_{nj} } - \sum_{\color{red}{h}=1}^{\color{red}{i}} \frac{ m_{hj} a_{hj}}{\left( 1+ \sum_{n=1}^{h} m_{nj} \right) \left( 2+\sum_{n=1}^{h} m_{nj} \right) }. $$
In this case, it seems that the approach suggested by @Kuifje for $c_{ij}$, cannot be applied for the $c^\prime_{ij}$ as the denominator of the second part of $c^\prime_{ij}$ is non-linear.
I have a couple of constraints over $M$ but all the constraints are linear.
Update 2: I think simulated annealing approach can be used for the cost function of $c^\prime_{ij}$.
Yes, linear programming is definitely a good approach for this problem.
Your objective function is
$$ f(m_{ij})=c_{ij}\;m_{ij} = \frac{a_{ij}{m_{ij}}}{1+\sum_{n=1}^i m_{nj}}-\frac{b_{ij}{m_{hj}m_{ij}}}{1+\sum_{n=1}^h m_{nj}} $$
Let $$t_{ij}=\frac{1}{1+\sum_{n=1}^i m_{nj}}\quad \mbox{and}\quad \omega_{ij}=\frac{1}{1+\sum_{n=1}^h m_{nj}}$$
Rewrite your objective function: $$ f(m_{ij},t_{ij},\omega_{ij})=\frac{a_{ij}{m_{ij}}}{1+\sum_{n=1}^i m_{nj}}-\frac{b_{ij}{m_{hj}m_{ij}}}{1+\sum_{n=1}^h m_{nj}} = a_{ij}m_{ij}t_{ij} - b_{ij}m_{hj}\omega_{ij}m_{ij} $$
and use the following change of variables: $$ x_{ij} = t_{ij}m_{ij} \quad \mbox{and}\quad y_{ij} = \omega_{ij}m_{ij} $$
You end up with $$ f(x_{ij},y_{ij},t_{ij},\omega_{ij})=a_{ij}x_{ij}-b_{ij}y_{ij}m_{hj} $$ subject to: \begin{cases} t_{ij}+t_{ij}\sum_{n=1}^i m_{nj}=1\\ \omega_{ij}+\omega_{ij}\sum_{n=1}^h m_{nj}=1 \end{cases} i.e., \begin{cases} t_{ij}+\sum_{n=1}^i x_{nj}=1\\ \omega_{ij}+\sum_{n=1}^h y_{nj}=1 \end{cases}
As you pointed out, the last step is to linearize the product $ y_{ij}m_{hj} $, which is fairly easy given that $y_{ij}$ is bounded by $1$ and $m_{hj}$ is a binary variable: replace $y_{ij}m_{hj}$ by $z_{ij}^h$ in the objective function and add the following constraints: \begin{cases} z_{ij}^h\le m_{hj}\\ z_{ij}^h\le y_{ij}\\ z_{ij}^h\ge y_{ij}+m_{hj}-1 \\ z_{ij}^h\ge 0 \end{cases}