Summary of the problem
---- 42 PARAMETER data %%% data set with 12 items cost hours people Rome 1000 5 5 Venice 200 1 10 Torin 500 3 2 Genova 700 7 8 Rome2 1020 5 6 Venice2 220 1 10 Torin2 520 3 2 Genova2 720 7 4 Rome3 1050 5 5 Venice3 250 1 8 Torin3 550 3 8 Genova3 750 7 8
We want to optimize the following objectives:
- Maximize number of items selected
- Minimize total cost
- Minimize total hours
- Minimize total number of people needed
In addition we have a number of simple bound constraints:
- The total cost can not exceed $10000
- The total number of hours we can spend is limited to 100
- There are only 50 people available
Non-dominated solutions.
We say that when comparing two feasible solution points \(a\) and \(b\), \(a\) dominates \(b\) if:
In this part I will show the GAMS code, that implements this.
We need to specify the problem data:
I will implement the implications as big-M constraints. This means we want to have small values for our big-M's. To obtain these it is useful to form tight bounds on the objectives \(f_k\).
- All the objectives of \(a\) are better or equal to the objectives of \(b\) and
- for one objective \(k\) we have \(f_k(a)\) is strictly better than \(f_k(b)\).
This means \(a\) does not dominate \(b\) if \(f_k(b)\) is strictly better than \(f_k(a)\) in one objective \(k\). Or if all objectives are the same.
Idea
The idea is to generate MIP solutions that have a preference to go to the efficient frontier using some weights, and then add constraints that say: we should be better in at least one objective. I.e. suppose we have collected already solutions \((\bar{x}_{i,s}, \bar{f}_{k,s})\) in previous cycles \(s\), then we require:\[ \begin{align} & \delta_{k,s} = 1 \Rightarrow f_k \text{ is better than } \bar{f}_{k,s} \\ & \sum_k \delta_{k,s} \ge 1 \>\> \forall s\\ & \delta_{k,s}\in \{0,1\}\end{align}\] I.e. our new solution is not dominated by previous record solutions.
This approach is documented in [3].
Note: I think we are just handwaving a bit about the possibility two solutions with all the same objectives: these are also non-dominated. For our problem this turns out not to be an issue.
Note: I think we are just handwaving a bit about the possibility two solutions with all the same objectives: these are also non-dominated. For our problem this turns out not to be an issue.
Approach 2: MIP models with non-dominated constraints
In this part I will show the GAMS code, that implements this.
2.1 Organize the data
We need to specify the problem data:
sets
i 'items' /Rome,Venice,Torin,Genova,Rome2,Venice2, Torin2,Genova2,Rome3,Venice3,Torin3,Genova3 / k 'objectives' /cost,hours,people,items/ ; table data(i,k) cost hours people Rome 1000 5 5 Venice 200 1 10 Torin 500 3 2 Genova 700 7 8 Rome2 1020 5 6 Venice2 220 1 10 Torin2 520 3 2 Genova2 720 7 4 Rome3 1050 5 5 Venice3 250 1 8 Torin3 550 3 8 Genova3 750 7 8 ; data(i,'items') = 1; parameter UpperLimit(k) / cost 10000 hours 100 people 50 /; * if zero or not listed use INF upperLimit(k)$(UpperLimit(k)=0) = INF; parameter sign(k) 'sign: +1:min, -1:max' / (cost,hours,people) +1 items -1 /; |
2.2 Form Box
I will implement the implications as big-M constraints. This means we want to have small values for our big-M's. To obtain these it is useful to form tight bounds on the objectives \(f_k\).
*---------------------------------------------
* Basic Model *--------------------------------------------- parameter c(k) 'objective coefficient'; variable z 'total single obj'; positive variable f(k) 'obj values'; binary variable x(i) 'select item' equations scalobj 'scalarized objective' objs(k) ; scalobj.. z =e= sum(k, c(k)*f(k)); objs(k).. f(k) =e= sum(i, data(i,k)*x(i)); * some objs have an upper limit f.up(k) = upperLimit(k); option optcr=0,limrow=0,limcol=0,solprint=silent,solvelink=5; *--------------------------------------------- * find lower and upper limits of each objective *--------------------------------------------- model m0 'initial model to find box' /all/; parameter fbox(k,*) 'new bounds on objectives'; alias(k,obj); loop(obj, c(obj) = 1; solve m0 minimizing z using mip; fbox(obj,'min') = z.l; solve m0 maximizing z using mip; fbox(obj,'max') = z.l; c(obj) = 0; ); display upperLimit,fbox; * range = max - min fbox(k,'range') = fbox(k,'max')-fbox(k,'min'); * tighten objectives f.lo(k) = fbox(k,'min'); f.up(k) = fbox(k,'max'); |
The first part is just the linear model: we have a scalarizing objective function, and the linear constraint that calculates the objectives \(f_k\).
In the second part we perform a loop, where we both minimize and maximize each \(f_k\).
It is interesting to see the differences between the implicit bounds we calculated here and the original constraints:
---- 81 PARAMETER UpperLimit cost 10000.000, hours 100.000, people 50.000, items +INF ---- 81 PARAMETER fbox new bounds on objectives max cost 6810.000 hours 45.000 people 50.000 items 9.000
We first observe that all lower bounds stay at zero. Actually a solution with \(x_i=f_k=0\) is a feasible and even non-dominated solution. However, on the upper bounds we did a really good job. These upper bounds were tightened significantly.
2.3 MOIP algorithm
*---------------------------------------------
* MOIP algorithm *--------------------------------------------- sets maxsols /sol1*sol1000/ sols(maxsols) 'found solutions' ; binary variable delta(maxsols,k) 'count solutions that are better'; equation nondominated1(maxsols,k) nondominated2(maxsols,k) atleastone(maxsols) cut(maxsols) ; parameter solstore(maxsols,*) M(k) 'big-M values' ; M(k) = fbox(k,'range'); scalar tol 'tolerance: improve obj by this much' /1/; * initialize solstore(maxsols,i) = 0; nondominated1(sols,k)$(sign(k)=+1).. f(k) =l= solstore(sols,k)-tol + (M(k)+tol)*(1-delta(sols,k)); nondominated2(sols,k)$(sign(k)=-1).. f(k) =g= solstore(sols,k)+tol - (M(k)+tol)*(1-delta(sols,k)); atleastone(sols).. sum(k, delta(sols,k)) =g= 1; cut(sols).. sum(i, (2*solstore(sols,i)-1)*x(i)) =l= sum(i,solstore(sols,i)) - 1; * initialize to empty set sols(maxsols) = no; * set weights to 1 (adapt for min/max) c(k) = sign(k); scalar ok /1/; model m1 /all/; loop(maxsols$ok, solve m1 minimizing z using mip; * optimal or integer solution found? ok = 1$(m1.modelstat=1 or m1.modelstat=8); if(ok, sols(maxsols) = yes; solstore(maxsols,k) = round(f.l(k)); solstore(maxsols,i) = round(x.l(i)); ); ); display solstore; |
In this algorithm we look for up to 1000 solutions. We implemented the constraints that make sure \(f_k\) is not-dominated by earlier solutions, and we use an objective that guides the solver towards the efficient frontier. In addition, for no good reason, I also added a constraint that works on the \(x\) space: forbid previously found integer solutions.
Note that in each cycle the MIP problem becomes bigger:
- we add non-dominated constraints
- we add an integer cut
- we add binary variables \(\delta_{s,k}\) related to the non-dominated constraints.
2.4 Solution
(Solution 0 with all zeroes is not printed here).
References
- Best Combination Algorithm of Complex Data with Multiple Constraints, https://stackoverflow.com/questions/55514627/best-combination-algorithm-of-complex-data-with-multiple-contraints
- Generating Pareto optimal points (part 1), https://yetanothermathprogrammingconsultant.blogspot.com/2019/04/generating-pareto-optimal-points-part-1.html
- John Sylva, Alejandro Crema: "A method for finding the set of non-dominated vectors for multiple objective integer linear programs", EJOR 158 (2004), 46-55.
No comments:
Post a Comment