Can this optimization problem be solved? 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}$.
 A: 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}
