Augmented form
With additional variables, we can have a more organized and simple set of
formulas, which should make the implementation process less prone to mistakes.
Of course, the downside would be the increase in model size and computation.
However, additional computational cost should be minimal.
Note
If this turns out to take too much time to compute, expanded form can be
used for optimization.
Now we formulate the problem as the following. Maximize
\[f(\mathbf{x}) = \sum_{i = 1}^{n_{\max}} p_{\mathrm{com}}(i) c_{i}
+ \sum_{i = 1}^{n_{\max}} \sum_{j = 1}^{j_{\max}(i)} \sum_{\sigma = 1}^{\sigma_{\max}}
p_{\mathrm{res}}(\sigma) r_{i, j, \sigma}\]
where
\[\begin{split}\mathbf{x} := & (\dots, c_{i}, \dots, r_{i}, \dots, e_{i}, \dots, s_{i}, \dots, r_{i, j}, \dots, e_{i, j}, \dots, s_{i, j}, \dots, r_{i, j, \sigma}, \dots, \\
& \dots, h_{i}, \dots, h_{i, j, \sigma}, \dots, l_{i, j, \xi}, \dots, u_{i}, \dots, \mu_{i}, \dots, \nu_{i}, \dots, \overline{a}_{i}, \dots a_{i}, \dots, \overline{b}_{i}, \dots, b_{i}, \dots)\end{split}\]
with
\[\begin{split}1 \leq i &\leq n_{\max} \\
1 \leq j &\leq j_{\max}(i) \\
1 \leq \sigma &\leq \sigma_{\max} \\
1 \leq \xi &\leq \xi_{\max}\end{split}\]
and
\[\begin{split}c_{i}, r_{i}, e_{i}, s_{i}, r_{i, j}, e_{i, j}, s_{i, j}, r_{i, j, \sigma}, \overline{a}_{i}, a_{i}, \overline{b}_{i}, b_{i} & \in \mathbb{R}_{0}^{+} \\
h_{i, j, \sigma}, l_{i, j, \xi}, u_{i}, \mu_{i}, \nu_{i} & \in \mathbb{Z}_{2} \\
h_{i} & \in \mathbb{N}_{0}\end{split}\]
subject to the real variable defining constraints
\[\begin{split}\begin{align}
& \forall i, & r_{i} - \sum_{j} r_{i, j} & = 0 \\
& \forall i, & e_{i} - \sum_{j} e_{i, j} & = 0 \\
& \forall i, & s_{i} - \sum_{j} s_{i, j} & = 0 \\
& \forall i, j & r_{i, j} - (e_{i, j} + s_{i, j}) & = 0 \\
& \forall i, j & r_{i, j} - \sum_{\sigma} r_{i, j, \sigma} & = 0 \\
& \forall i, & \overline{a}_{i} - \left(
(1 + \omega) (c_{i} + r_{i}) + \beta h_{i} + \kappa u_{i}
\right) & = 0 \\
& \forall i, & a_{i} - \left(
(1 + \omega) (c_{i} + e_{i}) + \omega s_{i} + \beta h_{i} + (\kappa - \tau\rho) u_{i}
\right) & = 0 \\
& \forall i, & \overline{b}_{i} - N \mu_{i} & \leq 0 \\
& \forall i, & \overline{b}_{i} - \overline{a}_{i} + N \mu_{i} & \leq N \\
& \forall i, & \overline{a}_{i} - \overline{b}_{i} + N \mu_{i} & \leq N \\
& \forall i, & b_{i} - N \mu_{i} & \leq 0 \\
& \forall i, & b_{i} - a_{i} + N \mu_{i} & \leq N \\
& \forall i, & a_{i} - b_{i} + N \mu_{i} & \leq N,
\end{align}\end{split}\]
the counting variable defining constraint
\[\begin{align}
& \forall i, & h_{i} - \sum_{j} \sum_{\sigma} h_{i, j, \sigma} & = 0
\end{align}\]
and the constraints for defining the relations between the two above, namely
\[\begin{split}\begin{align}
& \forall i, j, & r_{i, j} - N \sum_{\sigma} h_{i, j, \sigma} & \leq 0 & \\
& \forall i, j, & \sum_{\sigma} h_{i, j, \sigma} & \leq 1 & \\
& \forall i, j, \sigma, & N h_{i, j, \sigma} - r_{i, j} & \leq N - \delta_{\sigma, 0} - \varepsilon (1 - \varphi_{\sigma, 0}) \\
& \forall i, j, \sigma, & N h_{i, j, \sigma} + r_{i, j} & \leq N + \delta_{\sigma, 1} - \varepsilon (1 - \varphi_{\sigma, 1}) \\
& \forall i, j, \sigma, & r_{i, j, \sigma} - N h_{i, j, \sigma} & \leq 0 &
\end{align}\end{split}\]
for counting [redacted],
\[\begin{split}\begin{align}
& \forall i, j, & e_{i, j} - N \sum_{\xi} l_{i, j, \xi} & \leq 0 & \\
& \forall i, j, & \sum_{\xi} l_{i, j, \xi} & \leq 1 & \\
& \forall i, j, \xi, & N l_{i, j, \xi} - e_{i, j} & \leq N - \zeta_{\xi, 0} - \varepsilon (1 - \psi_{\xi, 0}) \\
& \forall i, j, \xi, & N l_{i, j, \xi} + e_{i, j} & \leq N + \zeta_{\xi, 1} - \varepsilon (1 - \psi_{\xi, 1})
\end{align}\end{split}\]
for counting [redacted],
\[\begin{split}\begin{align}
L_{\mathrm{com}} - \dot{L}_{\mathrm{com}} & \leq 0.5 \\
\dot{L}_{\mathrm{com}} - L_{\mathrm{com}} & \leq 0.5 - \varepsilon
\end{align}\end{split}\]
for counting [redacted], and
\[\begin{split}\begin{align}
& & u_{1} & = 1 & \\
& \forall i \geq 2, & (c_{i} + r_{i}) - N u_{i} & \leq 0 & (\dagger) \\
& \forall i \geq 2, & u_{i} - N (c_{i} + r_{i}) & \leq 0 & (\dagger) \\
& \forall i, & r_{i} - N \mu_{i} & \leq 0 & \\
& \forall i, & \mu_{i} - N r_{i} & \leq 0 & \\
& \forall i, & c_{i} - N \nu_{i} & \leq 0 & \\
& \forall i, & \nu_{i} - N c_{i} & \leq 0. &
\end{align}\end{split}\]
for counting [redacted], together with [redacted] constraints
\[\begin{split}\begin{align}
& & \sum_{i} a_{i} & \leq A_{B} \cdot R & \\
& & \sum_{i} \mu_{i} & \leq \mu_{\max}(T) & \\
& & \sum_{i} h_{i} & \leq h_{\max}(T) & \\
& & \sum_{i} b_{i} & \leq A_{F_{\mathrm{res}}}(T) & \\
& \forall i, & \overline{a}_{i} & \leq \overline{A}_{i} & \\
& \forall i \geq 2, & \overline{a}_{i + 1} - \overline{a}_{i} & \leq 0 & (\dagger) \\
& \forall i, & \mu_{i} + \nu_{i} & \leq 1 & \\
& \forall i \geq 2, & \nu_{i + 1} - \nu_{i} & \leq 0 & (\dagger) \\
& & \mu_{1} + \nu_{2} & \leq 1 & (\dagger) \\
& \forall i, j, & s_{i, j} - \alpha_{i} e_{i, j} & \leq 0 & \\
& & L_{\mathrm{res}} + L_{\mathrm{com}} & \leq P, &
\end{align}\end{split}\]
and user input custom constraints
\[\begin{split}\begin{align}
& \forall \sigma, & \sum_{i} \sum_{j} -h_{i, j, \sigma} & \leq -\hat{h}_{\min}(\sigma) \\
& \forall \sigma, & \sum_{i} \sum_{j} h_{i, j, \sigma} & \leq \hat{h}_{\max}(\sigma) \\
& & \sum_{i} -\mu_{i} & \leq -\hat\mu_{\min} \\
& & \sum_{i} \mu_{i} & \leq \hat\mu_{\max} \\
& & \sum_{i} -\nu_{i} & \leq -\hat\nu_{\min} \\
& & \sum_{i} \nu_{i} & \leq \hat\nu_{\max},
\end{align}\end{split}\]
and finally a set of computational cost optimizing constraint
\[\begin{align}
& \forall i, j, & r_{i, j + 1} - r_{i, j} & \leq 0.
\end{align}\]