## 3 Problem Statement
### Solving the Proposed Model  
The optimal schedule of a set of constellation tasks, including both individual and collective tasks, is achieved by solving the problem formally defined below:  
\[F: \max_{\mathbf{x}_{s,j,t}}\ \eqref{eq:constraint}\] (16) s.t. ( 3a )-( 15a )  
The resulting problem is an MILP (Mixed-Integer Linear Programming) problem, which has its general form given by:  
\[\min\ \mathbf{c}^{T}\mathbf{x}+\mathbf{h}^{T}\mathbf{y} \tag{17a}\] \[\text{s.t.}\ \mathbf{H}_{\mathbf{k}}\mathbf{x}+\mathbf{G}_{k} \mathbf{y}\leq\mathbf{b}_{k},\ \forall k\in\{1,\dots,m\}\] (17b) \[\ \ \ \ \ \ \mathbf{x}\in\mathbb{Z}^{n},\mathbf{y}\in\mathbb{R}^{p} \tag{17c}\]  
where the coefficients of the vectors \(\mathbf{c}\) and \(\mathbf{h}\) are easily obtained from (16), \(\mathbf{H}\) and \(\mathbf{G}\) are matrices and \(\mathbf{b}\) is a vector which are all obtained from the aforementioned \(m\) constraints.  
Typically, MILP problems can be solved to optimality or near optimality (within a dual gap) by exact algorithms such as branch-and-bound [22] or approximated using heuristics. In general, MILP is an NP-Hard problem in the sense that any NP-Complete problem (e.g., the classic satisfiability problem) can be reduced to it in polynomial time. Any polynomial time algorithm for MILP would imply that P = NP [A], in other words, all decision problems would be solved in polynomial time.  
In our work, instead of solving the problem in a straightforward manner defined by Equation (17), we employ the \(\epsilon\)-constrained methodology proposed by [PERSON] [20], which is illustrated in Figure 2. The methodology, employed for cost function (2), with terms A and B, is used to find the Pareto-optimal frontier of the proposed problem. From the perspective of the constellation operator, it is useful to give information about the trade-off that can be achieved by defining the priority for constellation and individual satellite tasks.  
The \(\epsilon\)-constrained method is known to provide efficient solutions for multi-objective mixed-integer problems. It consists of solving the problem multiple times along with a slack variable \(\epsilon\), which allows for the fixing of A or B values without rendering the problem infeasible. First, we maximize the constellation tasks (objective function term A), which provides an upper bound to the set of constellation tasks \(\mathcal{C}\) (\(A_{ub}\)). Then, we solve problem (17), now maximizing the individual tasks (objective function term B), given that the constellation tasks (A) QoS remains near the previously solved upper bound (within a given threshold); this step gives the lower bound for the individual tasks (\(B_{lb}\)). This step is then repeated, but now first maximizing B, then maximizing A, which will give the upper bound for the individual tasks (\(B_{ub}\)), and the lower bound of the constellation tasks (\(A_{lb}\)). These bounds are then considered in a next step as being the maximum and minimum range that both terms (A and B) can reach inside the Pareto-optimal frontier.  
In the next step, also illustrated in Figure 2, an iteration, guided by index \(j\), produces a Pareto-frontier with resolution of \(i+1\) points. The iterative process (over \(k\)) is responsible for calculating the Pareto points within the previously defined range \([A_{lb},A_{ub}]\), which relies on maximizing B in problem (17), once again ensuring that A remains within a desirable threshold, in an intermediate value between its maximum and minimum value:  
Figure 2: Representative flowchart of the \(\epsilon\)-constrained solution procedure.  
\[A=A_{lb}+k(A_{ub}-A_{lb})/i+\epsilon, \tag{18}\]  
whereby the obtained results are saved and used later to plot the Pareto-frontier.