Hybrid MPC with Bilinear Objective

[Behrooz]

Hi, 

I have a hybrid-MPC problem with a bilinear objective and some logic constraints implemented by "implies()" function. I use BMIBNB with fmincon as the upper solver and glpk as both the lower solver and LP solver for this problem, but I am not sure if this is the best way to handle this type of problems.

Some parts of my code is copied below. There are many more constraints in the original code, but since they are inherently the same as the constraints below I have not copied them here. As you can see in the code, the system dynamics is hybrid. Also different control laws are required for the two modes of the system. Finally, the cost function is bilinear in the output and input of the system (i.e. dTg(k) is the input and dx(3,k) is the output of the system).

I was wondering if there is a better way (in terms of calculation speed) for solving this problem (e.g. choosing a different solver).

Also, I was wondering if I can use TOMLAB/PENBMI as the solver rather than BMIBNB. 

Any help will be greatly appreciated.

%% Define the optimization variables

dTg = sdpvar(H,1);              

dx = sdpvar(nx, H);             

dBeta = sdpvar(H-1,1);          

Vend = sdpvar(1);

% Hybrid State dynamics
const = [];

const = [const, dx(:,1)==dx0];

for k = 1:H-1

    Model2 = [dx(:,k+1)==( A*dx(:,k)+B*dTg(k)+Bv*dV(k)+Bb*dBeta(k) ) ];

    Model3 = [dx(:,k+1)==( A130*dx(:,k)+B130*dTg(k)+Bv130*dV(k)+Bb130*dBeta(k) ) ];

    const = [const, implies(V(k)<=Vend-1e-3, Model2), implies(V(k)>=Vend, Model3)];

end

dxMax = [0.01; 0.01; 0.1];

const = [const, -repmat(dxMax,1,H)<=dx<=repmat(dxMax,1,H)];

for k = 1:H-1

    const = [const, implies( V(k)<=Vend, -1e-3<=dBeta(k)<=1e-3) ];

end

% Cost function

costfn = -sum( dx(3,1:H-1).* dTg(1:H-1) ) 

% Solve the problemoptions = sdpsettings('bmibnb.lowersolver','glpk','bmibnb.uppersolver','fmincon');

options = sdpsettings(options,'solver','bmibnb','bmibnb.lpsolver','glpk');

options = sdpsettings(options,'verbose',2,'solver','bmibnb');

options = sdpsettings(options,'debug',1);

out = solvesdp(const, costfn, options)

[Johan Löfberg]

Using penbmi instead of bmibnb makes no sense. Penbmi is a local solver for semidefinite programs.

Instead, what you perhaps is asking is

1. Is there a better local solver than fmincon to use in bmibnb

or

2. Is there a better global solver than bmibnb

answers
1. You can try ipopt, snopt or knitro.

2. YALMIP supports the global solvers baron and scip. Both most likely much better performing than bmibnb

and I hope you understand that you have a very very (insert more very) hard problem. Hybrid MPC is very hard. Nonconvex MPC is very hard. Hybrid MPC with nonconvex objective...

Finally, your hybrid model is not well posed. The solver can pick V(k) = Vend-1e-5, and then it is free to disregard both Model2 and Model3.

You are much better of by introducing an exclusive or construct

binvar d2 d3
[implies(d2,Model2), implies(d3,Model3), implies(d2, V(k)<=Vend), implies(d3,V(k)>=Vend), d2+d3 == 1]

[Johan Löfberg]

Have you skipped some code? V is not defined and is never really used in the model, it is only used to control other variables but is never really constrained it self?

[Behrooz]

Thank you very much for your fast response.

Yes, I have skipped some code. V is a known deterministic input to the system and its resolution ( V(k+1) - V(k) ) is known to be greater than 1e-3. 

Vend is an optimization variable that is chosen by the solver.

I avoided using binary variables since I thought that it will increase the complexity of the problem and make it more intensive.

I'll try the other options you have suggested.

Best,

[Johan Löfberg]

implies will introduce binary variables, so you are not saving anything by using it (implies is a logical statement, hence binaries required)

However, the model you have requires no binary variables. Since V is fixed, the model is trivially a linear prediction model. In fact, the command implies assumes the first argument is an sdpvar or constraint, and the model you have where the first argument is a constant, will lead to utter confusion in implies. Indeed, the implies operator crashes

sdpvar x
H = implies(5 >= 2, x <= 1);
Error in implies (line 49)

Are you sure V is constant? It should crash then

Anyway, if V is given, you should simply setup the linear model accordingly

for k = 1:H-1

 if V(k)<=Vend-1e-3

    const = [const,dx(:,k+1)==( A*dx(:,k)+B*dTg(k)+Bv*dV(k)+Bb*dBeta(k) ) ];
 elseif V(k)>=Vend

    const =[const, dx(:,k+1)==( A130*dx(:,k)+B130*dTg(k)+Bv130*dV(k)+Bb130*dBeta(k) ) ];
 else
    error('Behrooz has not defined this case');
 end

end

etc.

[Behrooz]

V is a constant vector, but Vend is a 1-by-1 sdpvar. 

[Johan Löfberg]

Aha. Then you must work with the exclusive or model.

[Johan Löfberg]

BTW, the structure should be exploited to add cuts on the binary variables

For instance, if V(i) <= V(j) you know that d2(i) must be true if d2(j) is true. Hence, a double for-loop with

if V(i) <= V(j)
 const = [const, d2(i) >= d3(j)]

etc.

[Behrooz]

I don't get it. Are you pursuing an approach to replace my logic implications by a double for-loop implementation as above? If yes, why is it required?

[Johan Löfberg]

First, you should not use your original implies constructions, since they are ill-posed (the solver can make the prediction undefined by placing Vend cleverly as I said). Instead, use the explicit exclusive or construct by defining two binary variables for every k and force them to sum to 1.

After that, it is, as always in integer programming, very beneficial if one can add additional redundant constraints on binary variables, which reduce the size of the continuous relaxation. In your case, you know a lot about the binary variables, since you know V (those are the cuts I told you to add)

Notice, although you could mistakenly think that you have exponential complexity 2^(H-1) in your model, that is not the case. Basically, there are only H combinatorial cases possible. They are

1. Vend <= min(V). All dynamics are the same

2. min(V) < Vend <= secondsmallestof(V). One dynamics follows model2, the rest model 3

3. etc

By using the O(H^2) cuts, you are helping the integer solver to discover this by making the continuous relaxation tighter.

Example to illustrate. I use d and 1-d instead of two separate binary variables to represent exclusive or. same thing. Let us look at the projection onto the d-variables for the continuous relaxation. You will see that it occupies a significant part of the hyper-cube

n = 3;
d = binvar(n,1);
v = randn(n,1);
vend = sdpvar(1)
const = [-10 <= vend <= 10];
for i = 1:n
    const = [const,implies(d(i), v(i) <= vend-1e-3)];
    const = [const,implies(1-d(i), v(i) >= vend)];
end
clf
plot(const,d,[],500,sdpsettings('relax',1))

Now add the cuts. Significant parts cuts, we're close to the case that only four corners are kept. i.e., the convex hull of the integer feasible points

for i = 1:n
    for j = i:n
        if v(i) <= v(j)
            const = [const, d(i) <= d(j)];
        end
    end
end
clf
plot(const,d,[],500,sdpsettings('relax',1))

I would advice you to solve these problems with and without cuts on your model, in an example where you use a linear or convex quadratic cost instead, to see if it has any benefit on those models when you are using a strong solver (if you keep the bilinear cost it will be hard to say anything since those problems are pretty much intractable and I think it would be hard to play around efficiently)

BTW, you should install a better MILP/MIQP solver, glpk is pretty much the worst one available.

[Behrooz]

you are right, my original logic implications were ill-posed. I have changed them according to your suggestion. However, by adding the redundant constraints to make cuts, more time is required to solve the problem.

Regarding the solver, I have CPLEX and have tried it as the LPSOLVER of bmibnb. The weird thing is that sometimes CPLEX as the LPSOLVER reports that the problem is infeasible while GLPK is able to solve the problem, and vice versa!

I appreciate it if you let me what you mean by a "better" solver. 

Thanks, 

[Johan Löfberg]

With better, I meant something like cplex, gurobi, mosek (as you used glpk in bmibnb, which indicated that you only had that solver available)

As I said, when you evaluate the use of cuts, do it on a MIQP or MILP version of your problem. The nonconvex part makes it really hard to say anything, since the bmibnb solver is so utterly basic (it might be that the sub-MILPS/sub-MIQPs in bmibnb are solved faster, but bmibnb it self has problems with the cuts for some reason)

How many nodes does cplex require on a MILP/MIQP version with and without cuts? 

If possible, post your model or send it to me.