Hybrid MPC Problem with constraint

[nms]

Dear Johan

I am developing a hybrid MPC in Simulink but I am currently having some difficulties with adding a certain constraint, please see attached file line 55-59.

The constraint I want to add is:

constraints = [constraints, Pb_min <= eta_Pout*x{k}(1)*u{k}(1) - eta_Pe*u{k}(2) <= Pb_max].

To avoid the non-linear term x{k}(1)*u{k}(1) I linearized it and obtained the following constraint instead  

constraints = [constraints, Pb_min <= eta_Pout*K*u{k}(1) - eta_Pe*u{k}(2) <= Pb_max, x1_lb <= x{k}(1) <= x1_ub],

where K depends on in between which boundaries x{k}(1) is found.

When running the simulation it runs for a while and then gets stuck inside the optimizer object. Once it is stuck I cant do anything and the only way (that I know of) to stop is by force quitting Matlab, thus I cant really analyse what happens when calling the optimizer.

I would appreciate it if you could give me some advice on what I am doing wrong or how I can solve the issue.

Please let me know if you would like all the files needed to run the simulation and I will email them to you.

Kind regards,

Niklas

[Johan Löfberg]

have you tried another solver? what is displayed? is cplex started? are you sending reasonable data to the optimizer object or are you sending inf/nans

[nms]

Just to add to the previous message. The simulation runs without any problems before I tried to add this constraint and the full simulation is about 800 seconds.

When I use cplex the simualtion runs for 70 seconds and then gets stuck inside the optimizer object and therefore I cant display any data. If I run the simulation for 60 seconds the results look good and the reference trajectory is followed.

If I use mosek the simulation runs for 4 seconds before the optimizer states that the problem is infeasible and returns nans. I have attached a file with the displayed output. I inserted breakpoints and stepped through the simulation and the input to the optimizer object seems reasonable.

I tried fmincon for a simulation time of 100 seconds and it manages to complete the simulation, but it is extremely slow. The results are pretty much garbage as the reference trajectory is not even close to being followed.

[Johan Löfberg]

When mosek thinks the problem is infeasible, what does cplex say for the same data?

Doesn't cplex display anything before it stalls? i.e, is it cplex that stalls before even starting to display stuff, or do you say it is optimizer which stalls?

[nms]

I tried to use the same input, that makes the problem infeasible when using mosek, for cplex and it does not have any trouble at all and returns a solution.

It is when the line of code below is run that the stall occurs when using cplex.

[solutions,diagnostics] = Controller(currentx,currentr,oldPe);

Unfortunately I cant save the complete simulation log since I cant press any button in Simulink but this is what is logged right before it stalls. 

Node log . . . Best integer = -3.464341e+06 Node = 0 Best node = -3.511324e+06 Best integer = -3.474102e+06 Node = 0 Best node = -3.511324e+06 Best integer = -3.511324e+06 Node = 0 Best node = -3.511324e+06 Implied bound cuts applied: 5 Flow cuts applied: 1 Mixed integer rounding cuts applied: 7 Tried aggregator 2 times. MIQP Presolve eliminated 417 rows and 34 columns. MIQP Presolve modified 238 coefficients. Aggregator did 3 substitutions. Reduced MIQP has 278 rows, 103 columns, and 904 nonzeros. Reduced MIQP has 56 binaries, 9 generals, 0 SOSs, and 0 indicators. Reduced MIQP objective Q matrix has 9 nonzeros. Presolve time = 0.01 sec. (1.30 ticks) Probing fixed 7 vars, tightened 14 bounds. Probing time = 0.00 sec. (0.27 ticks) Cover probing fixed 0 vars, tightened 85 bounds. Tried aggregator 2 times. MIQP Presolve eliminated 62 rows and 7 columns. MIQP Presolve modified 247 coefficients. Aggregator did 2 substitutions. Reduced MIQP has 214 rows, 94 columns, and 668 nonzeros. Reduced MIQP has 56 binaries, 0 generals, 0 SOSs, and 0 indicators. Reduced MIQP objective Q matrix has 9 nonzeros. Presolve time = 0.00 sec. (0.77 ticks) Probing fixed 2 vars, tightened 17 bounds. Probing time = 0.00 sec. (0.16 ticks) Cover probing fixed 0 vars, tightened 34 bounds. Tightened 1 constraints. MIP emphasis: balance optimality and feasibility. MIP search method: dynamic search. Parallel mode: deterministic, using up to 4 threads. Node log . . . Best integer = -3.377968e+06 Node = 0 Best node = -3.418581e+06

In Simulink it says "Running" but it has stopped at a certain time instant and there is no way of cancelling the simulation at this point. Pressing stop, pause, Ctrl+C, etc. does not do anything. I simply have to force quit Matlab at this point.

Thank you for taking your time to answer my questions, I appreciate it.

[Johan Löfberg]

OK, so cplex stalls.

Save all data and post it. I'm interested in both the Mosek failure and the cplex failure. To save a case, put a break and save to file using

save debugfile Controller currentx currentr oldPe

[nms]

Okay so I discovered that I had been a bit careless when I was trying out the optimizer using cplex with the currentx, currentr and oldPe that gave the optimizer using mosek an infeasible solution. Since oldPe showed as 1.000 in the workspace I simply typed oldPe = 1, but actually oldPe = 9.999958294000642e-01. I find this a bit strange though since oldPe is defined as a binvar. If I use Pe = 1 for the optimizer using mosek it will not return an infeasible solution.

For the problem regarding the stall of cplex have saved the values and controller that causes the stall, please see the cplexdebug.mat.

[Johan Löfberg]

solvers work with finite precision, and have tolerances on "integer enough". 

Which cplex version, and OS?

[nms]

Im running it on macOS 10.12 Sierra, Cplex 12.7, Matlab2016b. 

I have added the following lines to "definesolvers.m" in order to use Cplex 12.7.

solver(i) = qpsolver;

solver(i).tag     = 'CPLEX';

solver(i).version = 'IBM';

solver(i).subversion = '12.7.0';

solver(i).checkfor= {'cplexqcp.m','cplexlink1270'};

solver(i).call    = 'call_cplexibm_qcmiqp';

solver(i).objective.quadratic.nonconvex = 1;

solver(i).constraint.integer = 1;

solver(i).constraint.binary = 1;

solver(i).constraint.sos2 = 1;

solver(i).constraint.semivar = 1;

solver(i).constraint.semiintvar = 1;

solver(i).supportsinitial = 1;

i = i+1;

 

solver(i) = qpsolver;

solver(i).tag     = 'CPLEX';

solver(i).version = 'IBM';

solver(i).subversion = '12.7.0';

solver(i).checkfor= {'cplexqcp.m','cplexlink1270'};

solver(i).call    = 'call_cplexibm_qcmiqp';

solver(i).objective.quadratic.convex = 1;

solver(i).constraint.inequalities.elementwise.quadratic.convex = 1;

solver(i).constraint.inequalities.secondordercone.linear = 1;

solver(i).constraint.integer = 1;

solver(i).constraint.binary = 1;

solver(i).constraint.sos2 = 1;

solver(i).constraint.semivar = 1;

solver(i).constraint.semiintvar = 1;

solver(i).supportsinitial = 1;

i = i+1;

[Johan Löfberg]

I can only guess cplex has a memory bug on mac os then, unless you have an extremely memory-constrained machine

Only option is to switch solver, or randomly tweak options in cplex to see if you can get past this bug

[nms]

Regarding the problem I had when using mosek; if I add a Rouding Function block to the signal oldPe it does not run into the problem of the problem being infeasible. However the simulation is extremely slow compared to when I use cplex.

Unfortunately I am unable activate the license key for Gurobi due to issues with no reverse DNS information and the student version of Xpress seems to only be available for Windows.

Do you have any other suggestions of solvers I could try (academia)?

Thanks again for the help.

[Johan Löfberg]

I would try to tweak options in cplex to see if anything can be done. Turn off/change various features until the crash doesn't occur etc. Or downgrade to previous version perhaps and see if the bug is in there too. Randomly change the model (other ordering on constraints / variables)
 to see if anything changes in case it is a very weird bug

[nms]

Ok I will try to see if it works.

Stupid question maybe, but in which file can I tweak the options for cplex?

[Johan Löfberg]

An obvious tweak to your model, which also simplifies it, is to get rid of the implies for the definition of u

onstraints = [constraints,u{k}(2) == a{k}(1)*Pe_off + a{k}(2)*Pe_low+a{k}(3)*Pe_opt+*a{k}(4)*Pe_max]

or simply assign it before you use it as

u{k}(2) = a{k}(1)*Pe_off + a{k}(2)*Pe_low+a{k}(3)*Pe_opt+*a{k}(4)*Pe_max

[Johan Löfberg]

Besides the fact the cplex has a well-known but neglected bug which causes it to crash when the full options structure is sent via recent versions of matlab thus forcing us to use a limited set, you control it via sdpsettings

E.g

https://www.ibm.com/developerworks/community/forums/html/topic?id=79772c01-f545-44ef-9543-257197f1d9e7&ps=100

https://www.ibm.com/developerworks/community/forums/html/topic?id=0ed51b80-5898-4ccc-84f7-ce8ed2dca697&ps=100

[nms]

I tried both of your suggestions of simplification. The first suggestion gives the same results. The second suggestion, of assigning u{k}(2), does not run into any problems at all and runs at high speed through the whole simulation. Thank you very much! Hopefully I will not run into the same issue if I make tweaks to my model or constraints.