Let's consider the partial derivative of $an^2$ with respect to the variable $n_{i_0}$.
\begin{align*}
\frac{\partial\left(an^2\right)}{\partial n_{i_0}}
=\frac{\partial}{\partial n_{i_0}}\sum_{i}\sum_{j}n_in_j(1-k_{ij})\sqrt{a_ia_j}\tag{1}
\end{align*}
Note: We should avoid using the same variable name for different purposes. Since $i$ is used as bounded index variable on the RHS, we should not use $n_i$ as variable for the partial derivative on the LHS.
In order to calculate the partial derivative of double sums we have to do some bookkeeping by explicitely considering the parts of the double sum where the indices $i,j$ are equal with $i_0$.
Using for convenience the shorthand $A_{i,j}:=(1-k_{i,j})\sqrt{a_ia_j}$, we obtain from (1)
\begin{align*}
\frac{\partial\left(an^2\right)}{\partial n_{i_0}}&=\frac{\partial}{\partial n_{i_0}}\sum_{i}\sum_{j}A_{i,j}n_in_j\\
&=\frac{\partial}{\partial n_{i_0}}\sum_{i\neq i_0}\sum_{j\neq j_0}A_{i,j}n_in_j
+\frac{\partial}{\partial n_{i_0}}\sum_{i= i_0}\sum_{j= i_0}A_{i,j}n_{i}n_{j}\\
&\qquad+\frac{\partial}{\partial n_{i_0}}\sum_{i= i_0}\sum_{j\neq i_0}A_{{i},j}n_in_j
+\frac{\partial}{\partial n_{i_0}}\sum_{i\neq i_0}\sum_{j= i_0}A_{i,j}n_{i}n_{j}\tag{1}\\
&=\frac{\partial}{\partial n_{i_0}}\sum_{i\neq i_0}\sum_{j\neq i_0}A_{i,j}n_in_j
+\frac{\partial}{\partial n_{i_0}}A_{i_0,i_0}n_{i_0}^2\\
&\qquad+\frac{\partial}{\partial n_{i_0}}\sum_{j\neq i_0}A_{{i_0},j}n_{i_0}n_j
+\frac{\partial}{\partial n_{i_0}}\sum_{i\neq i_0}A_{i,{i_0}}n_{i}n_{i_0}\tag{2}\\
&=0+2A_{i_0,i_0}n_{i_0}+\sum_{j\neq i_0}A_{{i_0},j}n_j+\sum_{i\neq i_0}A_{i,{i_0}}n_{i}\tag{3}\\
&=\sum_{j}\left(A_{i_0,j}+A_{j,i_0}\right)n_j\tag{4}\\
\end{align*}
Comment:
In (1) we split the double sum in portions according to $i_0$ equal to $i,j$ resp. not equal to $i,j$.
In (2) we do some simplification. In fact (1) is only written for demonstration and we could start the apportionment with (2) instead.
In (3) we do the partial differentiation and rearrange the terms in (4)
If we substitute for $A_{i,j}$ back we obtain finally
\begin{align*}
\sum_{j}&\left(A_{i_0,j}+A_{j,i_0}\right)n_j\\
&=\sum_{j}\left((1-k_{i_0,j})\sqrt{a_{i_0}a_j}+(1-k_{j,i_0})\sqrt{a_ja_{i_0}}\right)n_j\\
&=\sqrt{a_{i_0}}\sum_{j}\left(2-k_{i_0,j}-k_{j,i_0}\right)n_j\sqrt{a_j}\tag{5}
\end{align*}
Note: If $k_{i,j}=0$ for all $i,j$ we obtain from (5) in accordance with OPs example
\begin{align*}
2\sqrt{a_{i_0}}\sum_{j}n_j\sqrt{a_j}
\end{align*}
Best Answer
It can be simplified somewhat.
$\begin{array}\\ f(N) &=\sum_{i=1}^{N}\sum_{j=1}^i \frac{1}{j^3}\\ &=\sum_{j=1}^N\sum_{i=j}^{N} \frac{1}{j^3}\\ &=\sum_{j=1}^N\frac{1}{j^3}\sum_{i=j}^{N} 1\\ &=\sum_{j=1}^N\frac{1}{j^3}(N-j+1)\\ &=\sum_{j=1}^N\frac{1}{j^3}(N+1)-\sum_{j=1}^N\frac{1}{j^3}j\\ &=(N+1)\sum_{j=1}^N\frac{1}{j^3}-\sum_{j=1}^N\frac{1}{j^2}\\ &\approx (N+1)\zeta(3)-\zeta(2) \qquad\text{for large }N\\ \end{array} $