This is a nonconvex objective with some linear constraints. Of course, the easiest way to gain some insight into how to solve this is to perform some quick and dirty tests with a global QP solver designed for these types of models. As the constraints are linear, a solver like Cplex or Gurobi may be a good starting point (pun intended).
If you have access to a local NLP solver, it may not always be easy to find a good starting point. One possible approach is to use a multistart algorithm (some NLP solvers, like Baron and Knitro, have this built-in). This will not guarantee a global solution (and most likely, it doesn't give you one), but at least we can prevent really bad solutions.
Model
It is very difficult to say much in the abstract. So let's create an example to put our teeth in. Here we have a transportation problem with a similar concave objective as in the original problem: \[\begin{align}\max&\sum_{(i,j)\in S}\color{darkred}x^2_{i,j} \\ & \sum_j \color{darkred}x_{i,j} \le \color{darkblue}{\mathit{supply}}_i && \forall i \\ & \sum_i \color{darkred}x_{i,j} \ge \color{darkblue}{\mathit{demand}}_j && \forall j \\ & \color{darkred}x_{i,j}\ge 0 \end{align}\] Here \(S\) is some subset of all \(\color{darkred}x_{i,j}\).
For completeness: this model fits the constraints of the original post: \(Ax \le b\), as we can write our constraints as \[\begin{align}& \sum_j \color{darkred}x_{i,j} \le \color{darkblue}{\mathit{supply}}_i && \forall i \\ & \sum_i \left(-\color{darkred}x_{i,j}\right) \le -\color{darkblue}{\mathit{demand}}_j && \forall j \\ & -\color{darkred}x_{i,j}\le 0\end{align}\]
Data generation
Using random supply and demand data, I generate a model with 50 demand nodes and 50 supply nodes and 2500 arcs (i.e. 2500 variables). The non-linear objective is determined by the set \(S\), which is in this case the left-upper submatrix of \((i,j)\). I use the following three sizes: the first 3, 10, and 20 rows and columns. This leads to 9, 100, and 400 non-linear terms.
I used Conopt (an excellent large-scale local NLP solver) with a few trials with different starting points. I have no clue about constructing good initial points, so I just generate random starting points. (In the answers to [1] it is suggested that finding good initial points is very easy -- I totally fail to see how that would or even could work.) The global QP solver I used is Cplex.
Size 1: 9 non-linear terms
This is a very easy problem. I used 10 multi-start iterations.
---- 162 PARAMETER res results (size=1) status obj time conopt .trial1 .* localopt 9572.357 0.032 conopt .trial2 . localopt conopt .trial3 . localopt 4711.833 0.016 conopt .trial4 . localopt 8382.324 0.031 conopt .trial5 .* localopt 9600.243 0.016 conopt .trial6 . localopt 3029.130 0.015 conopt .trial7 . localopt 0.016 conopt .trial8 . localopt 8035.934 0.031 conopt .trial9 . localopt 0.016 conopt .trial10. localopt 9466.108 0.015 multistart.total . 9600.243 0.188 cplex . . globalopt 9600.243 0.437
Note that the blank entries mean that the objective was zero. The records marked with an * are an improvement. We find the global optimal solution with conopt in trial 5.
Cplex is doing an excellent job here:
Nodes Cuts/
Node Left Objective IInf Best Integer Best Bound ItCnt Gap
0 0 9939.7624 0 9939.7624 7
* 0+ 0 9600.2432 9939.7624 3.54%
Found incumbent of value 9600.243195 after 0.02 sec. (23.55 ticks)
0 2 9939.7624 0 9600.2432 9939.7624 7 3.54%
Elapsed time = 0.38 sec. (2003.49 ticks, tree = 0.02 MB)
Size 2: 100 non-linear terms
---- 162 PARAMETER res results (size=2) status obj time conopt .trial1 .* localopt 19413.954 0.032 conopt .trial2 .* localopt 20443.936 0.032 conopt .trial3 . localopt 19391.875 0.031 conopt .trial4 . localopt 18962.226 0.016 conopt .trial5 . localopt 11144.281 0.016 conopt .trial6 . localopt 19914.037 0.031 conopt .trial7 . localopt 13040.606 0.015 conopt .trial8 . localopt 16662.541 0.031 conopt .trial9 . localopt 7131.035 0.016 conopt .trial10. localopt 18923.471 0.031 conopt .trial11. localopt 18073.915 0.032 conopt .trial12. localopt 11372.310 0.031 conopt .trial13. localopt 18614.723 0.031 conopt .trial14. localopt 17491.594 0.031 conopt .trial15. localopt 4861.867 0.015 conopt .trial16. localopt 19401.684 0.032 conopt .trial17. localopt 17091.372 0.031 conopt .trial18. localopt 14789.694 0.016 conopt .trial19. localopt 19026.705 0.032 conopt .trial20. localopt 16756.225 0.031 conopt .trial21. localopt 17776.559 0.016 conopt .trial22. localopt 13542.305 0.015 conopt .trial23.* localopt 20648.277 0.016 conopt .trial24. localopt 19895.247 0.031 conopt .trial25. localopt 16830.222 0.016 multistart.total . 20648.277 0.627 cplex . . globalopt 22491.252 61.641
The multistart algorithm finds a reasonable solution, but it is not optimal. There is about a 10% gap. The global solver needs more time to prove optimality.
Size 3: 400 non-linear terms
---- 162 PARAMETER res results (size=3) status obj time conopt .trial1 .* localopt 39316.482 0.047 conopt .trial2 .* localopt 44301.762 0.046 conopt .trial6 .* localopt 45559.234 0.047 conopt .trial16 .* localopt 45847.332 0.031 conopt .trial23 .* localopt 45922.772 0.031 conopt .trial34 .* localopt 46460.975 0.047 conopt .trial133.* localopt 46627.377 0.046 conopt .trial179.* localopt 47391.968 0.046 conopt .trial210.* localopt 47703.414 0.031 conopt .trial297.* localopt 48191.543 0.046 multistart.total . 48191.543 43.254 cplex . . nonopt 51650.218 1000.078
Here we used 1,000 multistart iterations. We only list the improvements here in the trace; otherwise, it would have 1,000 records. The global solver was stopped after hitting the time limit of 1,000 seconds. The gap at that time was 1.5%.
Multistart algorithm
The simplest multistart algorithm is also the one I used here:
- Create a random starting point
- Solve with a fast local solver
- If stopping criteria are met (e.g., time limit, iteration limit): STOP
- Go to step 1.
This is usually very fast, but we can't say much about global optimality. There are some techniques to say something about the probability that a better solution exists (this can be used as a more meaningful stopping criterion).
GAMS is not very well suited for large loops around a solve statement. Every time the model will be generated. This is costly. The GAMS scenario solver removes this burden by executing the loop inside the solver instead of GAMS itself. Unfortunately, the GAMS scenario solver does not support a multistart algorithms. The developers just forgot about this.
Global solvers
Deterministic global solvers behave much like MIP solvers. In fact, they use a branch-and-bound strategy, just like a MIP solver. We also can see this in how the "bounds" (best found solution so far vs. best possible bound) converge over time:
- Most of the excitement happens at the beginning: the solver finds quickly very good solutions, and the best bound is decreasing rapidly.
- The tail end is always a long slug: the solver is trying to prove optimality, mainly by working on the best bound.
- Often (but not always), the global optimal solution is found relatively quicky. But then, we need a lot of time proving optimality.
- In practice, we may set a time limit or use a stopping rule on the gap.
- We could have passed on the best multistart solution. This is automatic for global solvers like Baron. However, I don't know how to do that with Cplex. This is not a great loss, as Cplex often finds very good solutions very quickly. With a better starting point, we basically would have shaved some time from the early part of the solution process, where Cplex is already quite effective.
Data set
---- 32 PARAMETER supply supply1 17.175, supply2 84.327, supply3 55.038, supply4 30.114, supply5 29.221, supply6 22.405 supply7 34.983, supply8 85.627, supply9 6.711, supply10 50.021, supply11 99.812, supply12 57.873 supply13 99.113, supply14 76.225, supply15 13.069, supply16 63.972, supply17 15.952, supply18 25.008 supply19 66.893, supply20 43.536, supply21 35.970, supply22 35.144, supply23 13.149, supply24 15.010 supply25 58.911, supply26 83.089, supply27 23.082, supply28 66.573, supply29 77.586, supply30 30.366 supply31 11.049, supply32 50.238, supply33 16.017, supply34 87.246, supply35 26.511, supply36 28.581 supply37 59.396, supply38 72.272, supply39 62.825, supply40 46.380, supply41 41.331, supply42 11.770 supply43 31.421, supply44 4.655, supply45 33.855, supply46 18.210, supply47 64.573, supply48 56.075 supply49 76.996, supply50 29.781 ---- 32 PARAMETER demand demand1 69.574, demand2 79.046, demand3 66.208, demand4 31.850, demand5 12.106, demand6 13.715 demand7 67.589, demand8 57.995, demand9 6.616, demand10 82.700, demand11 10.740, demand12 21.030 demand13 56.027, demand14 78.485, demand15 21.276, demand16 6.878, demand17 61.977, demand18 65.587 demand19 42.400, demand20 39.335, demand21 27.767, demand22 28.106, demand23 16.514, demand24 96.809 demand25 41.458, demand26 81.804, demand27 33.467, demand28 16.012, demand29 78.351, demand30 10.387 demand31 23.665, demand32 3.970, demand33 30.425, demand34 53.449, demand35 18.592, demand36 20.881 demand37 36.528, demand38 35.154, demand39 35.672, demand40 99.861, demand41 102.824, demand42 40.454 demand43 40.753, demand44 80.662, demand45 43.132, demand46 94.773, demand47 15.422, demand48 77.012 demand49 9.006, demand50 61.094 ---- 32 PARAMETER totsup = 2245.137 total supply PARAMETER totdem = 2245.137 total demand ---- 32 SET subij subset to construct objective demand1 demand2 demand3 supply1 YES YES YES supply2 YES YES YES supply3 YES YES YES
Solution
---- 90 VARIABLE z.L = 9600.243 objective ---- 90 VARIABLE x.L flow from supply to demand nodes demand1 demand2 demand3 demand4 demand5 demand6 demand7 demand8 demand9 demand10 supply1 17.175 supply2 5.281 79.046 supply3 55.038 supply5 12.106 supply8 13.715 supply14 76.225 supply15 6.616 supply16 11.282 supply18 2.664 supply20 25.248 supply22 35.144 supply23 13.149 supply33 0.579 supply34 8.677 43.415 supply37 43.157 supply38 5.619 2.251 6.475 supply43 4.351 supply49 18.185 supply50 2.002 + demand11 demand12 demand13 demand14 demand15 demand16 demand17 demand18 demand19 demand20 supply4 16.730 supply8 37.154 11.917 supply12 19.299 supply19 24.911 supply24 6.005 supply28 50.059 supply30 21.276 9.090 supply32 48.857 supply36 25.111 supply38 39.335 supply39 1.731 supply41 41.331 supply47 33.310 supply50 10.740 6.878 + demand21 demand22 demand23 demand24 demand25 demand26 demand27 demand28 demand29 demand30 supply6 22.405 supply8 22.841 supply12 37.997 supply13 75.351 supply15 6.453 supply16 27.767 supply17 15.952 supply21 35.970 supply25 18.159 supply26 12.154 supply28 16.514 supply29 16.012 supply32 1.382 supply37 11.062 5.176 supply48 41.458 supply49 58.812 supply50 5.210 + demand31 demand32 demand33 demand34 demand35 demand36 demand37 demand38 demand39 demand40 supply5 1.693 supply10 23.665 supply11 35.161 11.028 supply13 23.763 supply20 18.287 supply26 70.935 supply34 35.154 supply36 3.470 supply38 18.592 supply40 9.852 36.528 supply45 3.430 30.425 supply47 31.263 supply50 0.540 4.410 + demand41 demand42 demand43 demand44 demand45 demand46 demand47 demand48 demand49 demand50 supply4 13.384 supply5 15.422 supply7 34.983 supply9 6.711 supply10 26.356 supply11 53.622 supply12 0.578 supply16 24.922 supply18 22.344 supply19 41.982 supply24 9.006 supply25 40.753 supply27 23.082 supply29 61.574 supply31 11.049 supply33 15.438 supply35 0.684 25.828 supply39 61.094 supply42 11.770 supply43 27.070 supply44 4.655 supply46 18.210 supply48 14.617
GAMS model
$onText
Nonconvex problem with linear inequality constraints
$offText
* three data sets
$set size 1
*$set size 2
*$set size 3
table sizedata(*,*) 'sizing'
NLsize QuadTerms MultistartIts CompactResults
size1 3 9 10 0
size2 10 100 25 0
size3 20 400 1000 1
* NLsize: size of left-upper (i,j) subset
* QuadTerms: NLsize^2
* MultistartIts: iterations for multistart algorithm
* CompactResults: 1=only list improvements
*---------------------------------------------------------------
* generate data for balanced transportation problem
*---------------------------------------------------------------
set
i 'supply nodes' /supply1*supply50/
j 'demand nodes' /demand1*demand50/
;
parameters
supply(i)
demand(j)
totsup 'total supply'
totdem 'total demand'
;
supply(i) = uniform(0,100);
demand(j) = uniform(0,100);
* make sure transportation problem is balanced
totsup = sum(i,supply(i));
totdem = sum(j,demand(j));
if (totsup>totdem,
demand(j) = demand(j) + (totsup-totdem)/card(j);
else
supply(i) = supply(i) + (totdem-totsup)/card(i);
);
totsup = sum(i,supply(i));
totdem = sum(j,demand(j));
display supply,demand,totsup,totdem;
*---------------------------------------------------------------
* a small subset has nonlinear terms in the objective
*---------------------------------------------------------------
scalar
NL 'square root of NL terms in obj'
;
NL = sizeData('size%size%','NLsize');
set subij(i,j) 'sub matrix';
subij(i,j) = ord(i)<=NL and ord(j)<=NL;
scalar nn 'number of non-linear terms';
nn = card(subij);
display nn;
*---------------------------------------------------------------
* concave nlp model
*---------------------------------------------------------------
Variables
z 'objective'
;
positive variables
x(i,j) 'flow from supply to demand nodes'
;
equations
obj
esupply
edemand
;
obj.. z =e= sum(subij,sqr(x(subij)));
esupply(i).. sum(j, x(i,j)) =l= supply(i);
edemand(j).. -sum(i, x(i,j)) =l= -demand(j);
model m /all/;
*---------------------------------------------------------------
* simple multistart algorithm
*---------------------------------------------------------------
parameter startx(i,j);
option qcp = conopt;
Sets
trialss 'superset' /trial1*trial10000/
trial(trialss) 'actual number of trials'
;
trial(trialss) = ord(trialss) <= sizeData('size%size%','MultistartIts');
* faster iterations
m.solprint = %solPrint.silent%;
m.solvelink = %solveLink.loadLibrary%;
parameter
res(*,*,*,*) 'results (size=%size%)'
best /-INF/
time /0/
improv 'improvement'
opt 'optimal'
print 'print this record'
printflag '1=all,0=only improvements'
;
option res:3:3:1;
singleton set imprv(*) 'improvement in trace';
printflag = 1-sizeData('size%size%','CompactResults');
acronym localopt,nonopt,globalopt;
loop(trial,
x.l(i,j) = uniform(1,100);
solve m maximizing z using qcp;
time = time + m.resusd;
opt = m.solvestat=%modelStat.Optimal% or m.solvestat=%modelStat.locallyOptimal%;
improv = opt and z.l>best;
print = printflag or improv;
best$improv = z.l;
if (print,
if (improv,
imprv('*') = yes;
else
imprv(' ') = yes;
);
res('conopt',trial,imprv,'status') = nonopt;
res('conopt',trial,imprv,'status')$opt = localopt;
res('conopt',trial,imprv,'obj') = z.l;
res('conopt',trial,imprv,'time') = m.resusd;
);
);
res('multistart','total',' ','obj') = best;
res('multistart','total',' ','time') = time;
*---------------------------------------------------------------
* global qp solver
*---------------------------------------------------------------
option qcp=cplex, reslim=1000;
m.optfile=1;
solve m maximizing z using qcp;
res('cplex',' ',' ','status') = nonopt;
res('cplex',' ',' ','status')$(m.solvestat=%modelStat.Optimal%)=globalopt;
res('cplex',' ',' ','obj') = z.l;
res('cplex',' ',' ','time') = m.resusd;
display res;
display z.l,x.l;
$onecho > cplex.opt
OptimalityTarget 3
$offecho
Conclusions
- I don't know how to generate a good starting point.
- A single run with a local NLP solver is not very useful: there is not much we can say about a solution found this way. It could be very bad or very good.
- The multistart approach is easy to implement, and finds reasonable solutions quickly. In general, it will not find global solutions.
- A global QP solver delivers solutions that are proven globally optimal, but it is expensive for larger problems. Often, we can stop a bit early, as much of the time is spent on proving optimality.
References
- Efficient Approaches for Separable Concave Quadratic Minimization over Linear Constraints?, https://or.stackexchange.com/questions/13081/efficient-approaches-for-separable-concave-quadratic-minimization-over-linear-co
No comments:
Post a Comment