You could rewrite the two last constraints as $$ f(i_1,i_2)-g(i_1,i_2) = e_{i_1i_2} $$ where $e_{i_1i_2}\in \mathbb{R}^+$ is a continuous variable that equals the difference between the left hand term and right hand term (and so $f(i_1,i_2)$ is your left hand term, and $g(i_1,i_2)$ your right hand term), and then minimize this variable in the objective function (eventually multiplied by a coefficient).
No, if you are working with linear programming, you cannot use the $\{<,>\}$ signs.
I think you need to work on a rolling horizon, that is, write your model for the the first month, and after week $1$, re-solve it. After week $2$, solve it again, and so forth. This of course implies that when you solve your problem, only the first week is definitely set. The other weeks are temporary and may change if balancing is required.
Although your model looks correct, have you considered other formulations such as a Dantzig-Wolfe decomposition? They fit your problem quite well :
Let $\mathcal{S}$ be the set of feasible shifts (e.g., over a week, or why not a month if it is relevant), and let $y_s$ be a binary variable that takes value $1$ if and only if shift $s\in \mathcal{S}$ is selected. Your problem then requires that :
- shifts are satisfied each day :
$$ \sum_{s \in \mathcal{S} \mid j \in s, t\in s} y_s = n_{j} \quad \forall j \in \{\mbox{morning},\mbox{afternoon},\mbox{night}\}, \quad \forall t=1,\dots,7 $$ (this means that for each day, and each shift type, you only sum on the shifts that have a working period on that day, for that shift type)
- you cannot have more than $\mid N \mid$ shifts ($1$ shift per nurse per week) : $$ \sum_{s \in \mathcal{S} } y_s \le|N| $$
To balance the shifts, you could define a coefficient $c_{s_is_j}$ for each pair of shifts that is more or less important if the shifts are pairwise balanced (e.g., $c_{s_is_j}=0$ if the shifts are balanced, $1$ if they are moderately unbalanced and $2$ if they are unbalanced). And so you could minimize $$ \sum_{s_i \in \mathcal{S}}\sum_{s_j \in \mathcal{S}} c_{s_is_j} \delta_{s_is_j} $$ where $\delta_{s_is_j}$ is another binary variable that takes value $1$ if and only if shifts $s_i, s_j$ are both selected (so you would have to add $y_{s_i}+y_{s_j} \le \delta_{s_is_j} +1$ to the model).
Of course this implies that you have generated your feasible shifts before hand (this is not too difficult). Alternatively, you could generate them dynamically if you are familiar with column generation.
