Not exactly "0" and "1" when solving MILP

[David Tena]

Hi! I am trying to solve a MILP using mosek. The point is that some of the sdpvar have logical operations on binvars such as "not" or "and", so exactly "0" and "1" solutions on binary variables are needed. Sometimes I get things like 8.4743e-17 instead of 0, and this is a problem because for matlab, not(0)=1 but not(8.4743e-17)=0. Is there something I can do? I have tried to use "round" command but I get an error. Thanks!

[Johan Löfberg]

Which solver?

10^-17 is simply numerical noise which makes absolutely no difference in practice. Just apply round (after you have extracted the solution)

[David Tena]

The solver is mosek. The problem is that there is relation between diferent sdpvars, as you can see in the next piece of code:

u=binvar(1,n-1);
x=sdpvar(1,n);

for k=1:n-1
    x(k+1)=0.85*x(k)+d(k)-e(k)*u(k); %State evolution
    if k>1
        c(k)=p(k)*u(k)+c_arranque*p(k)*and(u(k),not(u(k-1))); %Fails if u is 8.4743e-17 instead of 0
    else
        c(k)=p(k)*u(k);
    end
end

[David Tena]

Sorry, I am trying to minimize c

[David Tena]

I have discovered that if I use lpsolve instead of mosek, solutions of binvar have no numerical noise, they are always exacts 0 or 1, so my problem is solved with this. It seems to be a mosek problem. I don't know if it can be solved.

[G]

@David

how did you use lpsolve via Yalmip ?

I am pretty sure to have installed it correctly because if I try the examples provided it works, but I cannot set it as solver in Yalmip.

It tells me "solver no applicable lpsolve"

[Johan Löfberg]

What do you mean with fails?

If Mosek obtains a solution u(k)=8e-17, the expression and(u(k),not(u(k-1)) is surely 0. You can see that if you change your code to

temp(k) = and(u(k),not(u(k-1));
c(k)=p(k)*u(k)+c_arranque*p(k)*temp(k));

and then look at the value of temp. I have very hard to belive that anything else could arise, since the model for c=and(a,b) is [a>=c,b>=c,c >= a+b-1], so if any of those arguments are 1e-17, the solver would never let c be 0.

When you simulate it, simply round the obtained solution to get rid of the silly nuissance bits. Why Mosek leaves that is another issue

[Johan Löfberg]

So lpsolve rounds the solution. Don't let this silly thing convince you to switch from a state-of-the-art solver to a much weaker solver.

[Johan Löfberg]

Not applicable means you are trying to solve something else than a MILP

[G]

a = sdpvar(1,1);
options = sdpsettings('solver','lpsolve');
optimize(a>=1, 3*a, options)

[G]

a=intvar(1,1) of course, misprint in my last message

[Johan Löfberg]

Please show me the exact display when you ran that code

[David Tena]

El jueves, 5 de marzo de 2015, 18:16:34 (UTC+1), Johan Löfberg escribió:

What do you mean with fails?

If Mosek obtains a solution u(k)=8e-17, the expression and(u(k),not(u(k-1)) is surely 0.

That's what happens to me:

>> u=8e-17

u =

   8.0000e-17

>> and(u,1)

ans =

     1

>> 

[Johan Löfberg]

Well that is not the AND variable used by YALMIP. You are simply evaluating the AND operator on doubles. What I am saying is that internally, when the model is used by mosek, the variable modelling the AND value will have the value 0. Mosek never explicitly performs AND(uk,uk-1) if you think so. YALMIP creates a set of inequalities and variables such that the AND expression is represented (see last part of post)

I realized my example code will not work as intended due to the way AND is modelled vs evaluated a-posteriori. You have to introduce an auxilliary variable which simply copies the value of the AND operator, and look at that value instead (the value from AND is evaluated a-posteriori by simply evaluating the function on the values of the arguments, i.e, the way you do it)

temp1 = binvar(1,n);
temp2 = binvar(1,n);

...
temp1(k) = and(u(k),not(u(k-1));
c(k)=p(k)*u(k)+c_arranque*p(k)*temp1(k));
% Add to you list of constraints
Constraints = [Constraints, temp2(k) == temp1(k)]


% solve the problem

%and now look at the and operator value, i.e, this is the value the solver actually worked with
% when modelling and
value(temp2(k))

% This is not necessarily same as
and(value(u(k)),value(u(k-1))

since any numerical noise on the variables can destroy the logic



Alternatively said, if you simply solve an optimization problem

binvar a b
optimize(false(and(a,b)))

YALMIP writes this as

optimize([c <= a, c <= b, c >= a + b -1, c == 0])

A solver might give the solution a=1e-17, b = 1e-12 and c = 0. (numerical noise, all constraints satisfied). If you now evaluate the AND on the computed solution

and(value(a),value(b))

same thing happens if you evaluate the and operator (as YALMIP then simply evaluates the operator on the values of the arguments, it does not return the value of the variables used in the optimization

value(and(a,b))

and in both cases you will get 1. However, the solution computed by the solver is actually 0 which we can see if we introduce a dummy variable

binvar a b temp
optimize([false(and(a,b)),and(a,b) == temp])

the value of temp will be 0 (with possible noise)

You have to round any numerical noise if you suspect you have noise in binary variables, and you perform logic expression a-posteriori on the computed values

and(round(value(a)),round(value(b)))

[David Tena]

I have introduced the auxiliary variable as you say and now it works fine with Mosek. It is nice because until now it is the solver which has give me the best results. Thank you very much!

[Johan Löfberg]

It is not what I am saying you should do. You simply create a more complex model and you can just as well get these numerical effects on your new variables. To make your code robust, regardless of which solver you are using, you should clean the solution before you use it (using round)

[David Tena]

I thought I can only use ROUND command to get the final solutions when yalmip has finished its optimization, not in sdpvars ("c" is an sdpvar in my code). But you are right, if I use this code:

 temp(k)=round(and(u(k),not(u(k-1))));
 c(k)=p(k)*u(k)+c_arranque*p(k)*temp(k);

It works fine too

[Johan Löfberg]

It makes no sense to use round on a binary expression. It will simply add a redundant constraint and variable. If the model now happens to spit out a solution which is exactly integer, you are just seeing a random effect. You can probably obtain the same by other redundant changes to your model (reorder constraints, introduce dummy variables with no meaning etc epsilon changes to data)

[David Tena]

Ok, but if I use this:

temp(k)=and(round(value(u(k))),round(not(u(k-1))));
c(k)=p(k)*u(k)+c_arranque*p(k)*temp(k);

I get

Error using sdpvar/model (line 63)
Failed when trying to create a model for the "round" operator

Error in expandrecursive (line 19)
    [properties,F_graph,arguments,fcn] =
    model(variable,method,options,allExtStruct(ext_index),w);

Error in expandrecursive (line 206)
                            [F_expand,failure,cause] =
                            expandrecursive(recover(expressionvariables(i)),F_expand,extendedvariables,monomtable,variabletype,where,level+1,options,'exact',[],'none',allExtStruct,w);
                           
Error in expandmodel>expand (line 376)
            [F_expand,failure,cause] =
            expandrecursive(recover(variables(i)),F_expand,extendedvariables,monomtable,variabletype,where,level+1,options,method,[],'convex',allExtStruct,w);
           
Error in expandmodel (line 209)
    [F_expand,failure,cause] =
    expand(index_in_extended,variables,h,F_expand,extendedvariables,monomtable,variabletype,'objective',0,options,method,[],allExtStructs,w);
   
Error in compileinterfacedata (line 116)
    [F,failure,cause,operators] = expandmodel(F,h,options);

Error in solvesdp (line 251)
[interfacedata,recoverdata,solver,diagnostic,F,Fremoved,ForiginalQuadratics]
=
compileinterfacedata(F,[],logdetStruct,h,options,0,solving_parametric);

Error in optimize (line 31)
varargout{:} = solvesdp(varargin{:});

Error in PTAR_binvar_resthoras2 (line 50)
sol=optimize(rest,obj,opciones)

I think I am missing something...

[David Tena]

Sorry, I mean

temp(k)=and(round(u(k)),round(not(u(k-1))));
c(k)=p(k)*u(k)+c_arranque*p(k)*temp(k);

[Johan Löfberg]

I don't. Old version?

u = sdpvar(3,1);
k = 2;
temp(k)=and(round(value(u(k))),round(not(u(k-1))))
Nonlinear matrix variable 1x2 (full, real, binary, 1 variable, values in range [0,1])

However, this code makes no sense either. u is binary, so applying a round operator is entirely redundant and just leads to a model with more but redundant variables and constraints. YALMIP and the solver is working with all variables defined as integer/binary as precisely that. The fact that the solver returns a slightly non-binary solution is just a consequence of the fact that the solver works in finite precision and uses numerical optimization to solve the problem. If you must have an exactly binary variables and use a solver which doesn't round the end result for you, you have to apply round your self a-posteriori, i.e, once you have computed the solution and actually use it, round(value(u(k))

[Johan Löfberg]

Your model

binvar a b
optimize(round(a) == b)

Is simply extended to (c is an integer which is forced to 0 if a<.5, and 1 if a>.5 and undefined 0 or 1 when a=.5)

binvar a b
intvar c
optimize([c == b, a-0.5  <= c <= a + 0.5])

Hopefully you see that this just complicates the model, and adds absolutely no value since a only can take the values 0 or 1 anyway

[David Tena]

That's my entire code, sorry for the comments in spanish. I get the error I mentioned before. Don't you?

poolh=[30.01 25.05 20.43 13 12.4 16.5 26.11 43 40 48 44 39.6 38.4 34.51 28.63 28.6 34.2 37.08 37.59 48.8 51.48 43.1 38.55 37.33];
n=max(size(poolh)); %datos por día
MWsoplantes=0.5; %Dos soplantes de 250 kW según manual
precioh=poolh+48.371-0.5-0.07*(poolh+48.371); %MWh
p=MWsoplantes*precioh; %precio de tener en marcha las soplantes cada hora

%Constantes
e=1.5*ones(1,n); %Tasa máxima de eliminación de amonio con soplantes al máximo
d=2*ones(1,n); %Pendiente de incremento de amonio si no se airea. Función del flujo de amonio entrada
x_ini=9; %Amonio inicial

%Restricciones
hsin=4; %Número de horas máximas sin soplar
c_arranque=2; %Coeficiente de consumo eléctrico extra por arranque (% de potencia extra consumida)
xlimit=10; %Amonio máximo permitido

% Optimizar u
u=binvar(1,n-1);
x=sdpvar(1,n);
horas_sin=sdpvar(1,n);
temp=binvar(1,n-1);

for k=1:n-1
    x(k+1)=0.85*x(k)+d(k)-e(k)*u(k); %Evolución del estado
    if k>1
        %temp(k)=round(and(u(k),not(u(k-1))));
        temp(k)=and(round(u(k)),round(not(u(k-1))));
        c(k)=p(k)*u(k)+c_arranque*p(k)*temp(k);
    else
        %temp(k)=round(u(k));
        c(k)=p(k)*(1+c_arranque)*u(k);
    end
   
    if k>=hsin+1 %Bucle para contar horas seguidas sin soplar
        horas_sin(k)=0;
        for j=1:hsin+1
            horas_sin(k)=horas_sin(k)+u(k-j+1);%Variable "hsin" últimas horas (sin soplar)
        end
    end
end

rest=[x(1)==x_ini, x<=xlimit , x>=0, horas_sin>=1];
obj=ones(1,n-1)*c';
opciones=sdpsettings('solver','mosek','verbose',0);

[David Tena]

sol=optimize(rest,obj,opciones)

[Johan Löfberg]

Yes I do, when I actually call the solver

>> optimize(rest, obj)

Error using sdpvar/model (line 63)
Failed when trying to create a model for the "round"
operator

Interesting. However, luckily not an issue for you, since you shouldn't use the round operator

[David Tena]

I actually need the round operator because if u(k) or u(k-1) has numerical noise, the results of:

and(u(k),not(u(k-1)))

change.

For exemple, execute the next code. The solution of the variable "cost" change if you use the two different lines of code for temp(k) (first one is now commented). The right solution (the one who does what I need) is the one that I get with the second line (using round). If I don't use round, you can see in the figure that the first green column is not so tall.


clear all
clc
tic
%Datos costes, pool
poolh=[30.01 25.05 20.43 13 12.4 16.5 26.11 43 40 48 44 39.6 38.4 34.51 28.63 28.6 34.2 37.08 37.59 48.8 51.48 43.1 38.55 37.33];

n=max(size(poolh)); %datos por día
MWsoplantes=0.5; %Dos soplantes de 250 kW según manual
precioh=poolh+48.371-0.5-0.07*(poolh+48.371); %MWh
p=MWsoplantes*precioh; %precio de tener en marcha las soplantes cada hora

%Constantes
e=1.5*ones(1,n); %Tasa máxima de eliminación de amonio con soplantes al máximo
d=2*ones(1,n); %Pendiente de incremento de amonio si no se airea. Función del flujo de amonio entrada
x_ini=9; %Amonio inicial

%Restricciones
hsin=4; %Número de horas má

ximas sin soplar
c_arranque=0.2; %Coeficiente de consumo eléctrico extra por arranque (% de potencia extra consumida)
xlimit=10; %Amonio máximo permitido

% Optimizar u
u=binvar(1,n-1);
x=sdpvar(1,n);
horas_sin=sdpvar(1,n);
temp=binvar(1,n-1);

for k=1:n-1
    x(k+1)=0.85*x(k)+d(k)-e(k)*u(k); %Evolución del estado
    if k>1

        %temp(k)=and(u(k),not(u(k-1)));
        temp(k)=round(and(u(k),not(u(k-1))));

        c(k)=p(k)*u(k)+c_arranque*p(k)*temp(k);
    else
        %temp(k)=round(u(k));
        c(k)=p(k)*(1+c_arranque)*u(k);
    end
   
    if k>=hsin+1 %Bucle para contar horas seguidas sin soplar
        horas_sin(k)=0;
        for j=1:hsin+1
            horas_sin(k)=horas_sin(k)+u(k-j+1);%Variable "hsin" últimas horas (sin soplar)
        end
    end
end

rest=[x(1)==x_ini, x<=xlimit , x>=0, horas_sin>=1];
obj=ones(1,n-1)*c';
opciones=sdpsettings('solver','mosek','verbose

',0);

sol=optimize(rest,obj,opciones)
cs=double(c);
cost=ones(1,n-1)*cs'

toc
xs=double(x);
xlim=xlimit*ones(1,max(size(xs)));
us=round(double(u));

clf
stairs(p,'k')
hold on
plot(xs,'b')
hold on
plot(xlim,'b --')
hold on
stairs(cs,'g')
hold on
stairs(us,'r')
legend('Precio EE','Concentración amonio', 'Concentración máx','Coste','Soplantes')
grid

[Johan Löfberg]

But you are doing precisely the mistake I am talking about. When c is evaluated a-posteriori, it does not extract the value of the AND operator mosek worked with when YALMIP had sent a model of the AND operator, but evaluates the AND operator for the values of u returned by mosek. If those have noise, the c used while solving, and the c evaluated for particular noisy solution, might differ

Hence, either round the solution before doing anything

assign(u,round(value(u)));

and you will have the same regardless of the use of round.

Or, as a test to see what really is computed, evaluate everything from scratch

clear x
clear c
clear temp
x = x_ini;
u = round(value(u));

for k=1:n-1

    x(k+1)=0.85*x(k)+d(k)-e(k)*u(k); %Evolución del estado

    if k>1       

        temp(k)=(and(u(k),not(u(k-1))));

[Johan Löfberg]

Note, the use of the round operator will indeed "fix" the evaluation of the expression in c, but it has not in any sense fixed the values of u. Hence, if you evaluate from scratch (i.e, if you actually want to use the solution for something), you still have to use round on your solution. The use of round in the model is just a band-aid, it does not solve the underlying problem.

[David Tena]

Right. ¿What do you think about this solution? It seems to work well and no ROUND is used. I have treated temp as a constraint.

% Restricción de horas sin soplar. Incremento en el coste de arranque
clear all
clc
tic

%Datos costes, pool
poolh=[30.01 25.05 20.43 13 12.4 16.5 26.11 43 40 48 44 39.6 38.4 34.51 28.63 28.6 34.2 37.08 37.59 48.8 51.48 43.1 38.55 37.33];
n=max(size(poolh)); %datos por día
MWsoplantes=0.5; %Dos soplantes de 250 kW según manual
precioh=poolh+48.371-0.5-0.07*(poolh+48.371); %MWh
p=MWsoplantes*precioh; %precio de tener en marcha las soplantes cada hora

%Constantes
e=1.5*ones(1,n); %Tasa máxima de eliminación de amonio con soplantes al máximo
d=2*ones(1,n); %Pendiente de incremento de amonio si no se airea. Función del flujo de amonio entrada
x_ini=9; %Amonio inicial

%Restricciones
hsin=4; %Número de horas máximas sin soplar
c_arranque=0.2; %Coeficiente de consumo eléctrico extra por arranque (% de potencia extra consumida)
xlimit=10; %Amonio máximo permitido

% Optimizar u
u=binvar(1,n-1);
x=sdpvar(1,n);
horas_sin=sdpvar(1,n);
temp=binvar(1,n-1);

rest=[];

for k=1:n-1
    x(k+1)=0.85*x(k)+d(k)-e(k)*u(k); %Evolución del estado
    if k>1
        %temp(k)=and(u(k),not(u(k-1)));

        %temp(k)=round(and(u(k),not(u(k-1))));
        c(k)=p(k)*u(k)+c_arranque*p(k)*temp(k);
        rest=[rest, temp(k)==and(u(k),not(u(k-1)))];
    else
        c(k)=p(k)*(1+c_arranque)*u(k);
        rest=[rest, temp(k)==u(k)];

    end
   
    if k>=hsin+1 %Bucle para contar horas seguidas sin soplar
        horas_sin(k)=0;
        for j=1:hsin+1
            horas_sin(k)=horas_sin(k)+u(k-j+1);%Variable "hsin" últimas horas (sin soplar)
        end
    end
end

rest=[rest, x(1)==x_ini, x<=xlimit , x>=0, horas_sin>=1];

[Johan Löfberg]

Sure, but it is all band-aids to come around the issue that you can have numerical noise in solutions, and with the discontinuous behavior of logic expressions, one has to be careful when working with the solutions. You still have to clean u if you want to use it in practice.

[G]

What do you mean by 'display'

[Johan Löfberg]

Everything printed in the MATLAB window when you run the code

[G]

>> a = intvar(1,1)
Linear scalar (real, integer, 1 variable)
>> options = sdpsettings('verbose','1','solver','lpsolve');
>> optimize(a>=1, 3*a, options)
Warning: Solver not found (lpsolve)
ans =
    solvertime: 0
          info: 'Solver not found (lpsolve)'
       problem: -3
    yalmiptime: 0.0945
>>

[Johan Löfberg]

So? This is not what you said above, and the problem is precisely what the error message says

[Johan Löfberg]

http://users.isy.liu.se/johanl/yalmip/pmwiki.php?n=Solvers.LPSOLVE