Difficulty finding the real issue with a constraint

[Tom Cunningham]

Hello Dr. Lofberg,

I'm on my second project using YALMIP / Gurobi.  The first one went very well, but I'm having difficulty with my current project, which is writing a MILP problem to determine the optimum sequence of spacecraft operating modes at every timestep.  Currently a 100-minute simulation, approximately half an orbital period, one minute per time step.

I'm using for my decision variable 

mode = binvar(J,K);

Where J is the total number of time steps and K is the total number of possible modes.  It may be possible, at a given j value, for there to be multiple k's set to 1, that is, multiple modes active, but for now I'm constraining it so just one k be equal to 1. 

My battery constraints are setup to track battery state at every minute, and I do it step-wise.  So if the battery state is known for jj = 1, then for every jj from jj = 2 to 100, I use a loop...

for jj = 2:J

    % 2. Build batt_state_inst(jj), set it >= batt_min_charge

 

    batt_state_inst(jj) = batt_state_inst(jj-1) + mode(jj,1)*safe_powerrate(jj)  + mode(jj,2)*hires_powerrate(jj) + mode(jj,4)*downlink_powerrate(jj) ...

    + mode(jj,3)*imu_powerrate(jj) + mode(jj,5)*safe_nocharge(jj);

 

    F = [F, batt_state_inst(jj) >= batt_min_charge];

etc.

So the instantaneous batter state is built at time step jj, as the previous state plus the mode variable for that time step, which would result in the battery charge being increased or decreased.  I was hoping this would be a valid setup for YALMIP constraints, but I'm getting the error 

Error using lmi/horzcat (line 21)

One of the constraints evaluates to a FALSE LOGICAL variable

Error in reducedscaleversion (line 77)

    F = [F, batt_state_inst(jj)-batt_min_charge >= 0]; 

I've tried a few variations on the code but always get the same error.

Could it be that I need to use the implies command?  How would I set that up for multiple mode variables (all mode(jj,:) for a given jj), each with its own power usage?

Or perhaps there is another issue I'm not seeing?

Many thanks,

Tom

[Johan Löfberg]

It means batt_state_inst(jj)-batt_min_charge evaluates to a negative constant, which trivially cannot be >=0. Evaluate batt_state_inst(jj)-batt_min_charge in the code and you will see it displays a constant before it crashes

Perhaps J=K and thus you get a symmetric variables when you want a fully parameterized?

BTW, I hope mode is the only decision variable in the expression when you  generate the trajectory batt_state_inst. Otherwise it is horribly nonlinear. How is batt_state_inst initialized? This looks like code which typically is incorrectly coded as you cannot do state(1)=constant, state(2) = state(1) + decisionvariables with objects in matlab. The preferred pattern is statetrajectory = [statetrajetory nextstate]

[Tom Cunningham]

Thank you, I think you may have found the issue.

As for the decision variable mode(J,K), J is not equal to K, but to be sure I changed the definition to 

mode = binvar(J,K,'full');

just in case, but it did not change the error.   Yes, this is my only decision variable here.  

You correctly summarized that my code is trying to do this pattern:

>>  state(1)=constant, state(2) = state(1) + decisionvariables 

...and you say this can't or shouldn't be done in matlab. 

Could you please elaborate a little on the preferred pattern...  What does that pattern do differently?  How might that be applied to the batt_state expression?

Thanks again for your guidance!

[Johan Löfberg]

also note that a sparser and typically just as good or even better model is

batt_state_inst = sdpvar(J,1);
Model = [batt_state_inst (2:end) = batt_state_inst(1:end-1)  + mode*rates, batt_state_inst(1) == givenstate]

(where you stacked all the rate vectors in the same matrix to get the code working)

[Tom Cunningham]

Interesting...define another decision variable for batt_state, making use of the previously-defined mode binvar.

I'll give it a try.  Cheers!

[Johan Löfberg]

sdpvar x
y(1) = 1;
y(2) = y(1)  + x

simply cannot be done in matlab with new classes, as built-in has precedence in assignment, so the code is equivalent to 

sdpvar x

y(1) = 1;

y(2) = y(1)  + double(x)

which thus will evaluate to NaN as x doesn't have any value yet ( and would make sense even if x had a value)

hence, we do

y(1) = 1
y= [y y(1)+x]

instead, or variants of that

[Johan Löfberg]

it's the standard way to represent dynamical systems optimization problems

https://yalmip.github.io/example/standardmpc/