How can I use 'OR' constraints with fmincon / linprog: { x=0 } or { 0.6 <= x <= 1 }

Hello!

I have a matlab script that optimize the electricity production from a heating plant based on the current electricity price and available boilers. Its very close to: http://www.mathworks.se/products/optimization/description6.html

However, the boilers can only run between 60-100 % or be completely turned off. So the constraints should be like:

x=0

or

0.6 <= x =< 1

for every boiler, and I have 10 of them.

How can I put this into linprog or fmincon? I want to have a script that tells me when to turn the boilers on and off and also at what rate every boiler should go.

Thank you very much

Answer 1

[Matt J]

The only robust way is to loop over all 2^10 on/off combinations of the boilers. In each pass through the loop, i.e., for each combination, use fmincon/linprog to solve for the optimal output of the on-boilers. Shouldn't be too bad with only 10 unknowns max. per problem.

[jar]

Thank you very much for the quick respond!

The workaround you just suggested is exactly what I have done previously and its OK for small amount of calculations (takes around 5-6 seconds). But I want to add a storage tank to calculate what's best for a period of 24 hours and it scales very fast.

I'm trying to get this working:

http://www.mathworks.se/help/optim/ug/writing-constraints.html#btqvfel

but I cant understand how to do it with 10 variables :/

[Matt J]

The technique there only has a chance of working if the constrained region is connected which, in your case, it is not. And this still ignores the risks of the non-smoothness of c(x), mentioned at the link.

You could reformulate your problem in terms of expressions b(i)*x(i) where x(i) is boiler output and b(i) is a binary variable describing the state of the i-th boiler, (1=on, 0=off). If you have the Global Optimization Toolbox, you could then try applying GA to the reformulated problem. However, it's hard say whether this will do any better for you than exhaustive search.

[jar]

Ah ok, thank you!

I think I found another way of doing it without the GA toolbox, see below. It works fine for almost every case but not everytime, but I have no idea why :(

Answer 2

[Matt J]

If your objective is linear, another idea might be to discretize x(i) and use bintprog(). For example, using binary variables y(i,j), the following

 x(i)=0.6*y(i,1)+0.1*(y(i,2)+y(i,3)+y(i,4)+y(i,5));

       s.t.  4*y(i,1)>=y(i,2)+y(i,3)+y(i,4)+y(i,5),

allows x(i) to take on values [0,.6,.7,.8,.9,1.0]. You could use a bintprog solution to this approximate problem to first decide which boilers should be off/on. Once you know which biolers are on/off, you can then solve for the on-boiler x(i) more precisely and in continuous space using linprog.

[jar]

I will try that, thank you very much!

Answer 3

[jar]

Hi Matt! I think I found a way to do it with fmincon without using bintprog or GA.

As you know I want x(1) to be my boiler input, so that x(1) is equal to 0 or in between 0.4-1.

Let me introduce three help variables x(3), x(4) and x(5).

x(3) can vary between 0.4 and 1.

x(4) and x(5) can be equal to any digit between 0 and 1. I set these constraints in lb and ub.

So the following constraints will be set as nonlinear equalities constraints.

x(4)*x(5)=0 , this means either x(4) or x(5) has to be zero.

After that, I put the following in Aeq:

x(4)+x(5)-x(3)=0, this means that the non zero variable of either x(4) or x(5) has to be equal to x(3), (0.4 - 1).

And Finnaly:

x(1)=(x(3)+x(4)-x(5))/2 which will give me exactly what I'm aiming for. If x(4)=0, then x(1)=0 and if x(5)=0, then x(1)=x(3), in other words somewhere between 0.4 and 1.

However, this works fine for all values of power demand in my test program eccept for when the demand is equal to 5?? I have also introduced an electric heater that can run anywhere between 0 to 10 MW to cover up when the demand is < 6. But I have no clue why this isnt working when the power demand is equal to 5? It works fine otherwize, for example when the power demand is 4 or 6.

function test
close all,clear all,clc

%%Input

load=5 % [MW]

x1=15; % Maximal boiler output [MW]

x2=10; % Electric heater [MW]

% help variables for the boiler

x3=x1;

x4=x1; 

x5=x1;

%%Bounds

ub=ones(1,5); % Upper bounds

lb=zeros(1,5); % Lower bounds

lb(3)=6/x1; % Boiler have to be equal to 0 or in between 0.4 - 1.

Aeq=[0 0 x3 -x4 -x5; % x(4) or x(5) must be equal to x(3), see f. nonlcon

      -x1 0 x3/2 x4/2 -x5/2; % x(1) is 0 or equal to x(3)

      x1 x2 0 0 0]; % The boiler and/or the elecric heater = load

beq=[0
    0
    load];

%%Optimize!

x0 = [0 0 1 0 0]; 

options = optimset('Algorithm','active-set');

[xval fval]=fmincon(@myfun,x0,[],[],Aeq,beq,lb,ub,@nonlcon,options);

x=[x1 x2 x3 x4 x5];

xval.*x

end

%%Objective function

function [f] = myfun(x)

f=x(1)+10*x(2); % The electric heater is 10 times mor expensive to run

end

%%Nonlinear constrains

function [c ceq] = nonlcon(x)

c=[];

ceq(1)=x(4)*x(5); % Either x(4) or x(5) has to be equal to 0

end

[Matt J]

Hi jar,

It can't reliably work unfortunately. The issue was never about how to express the constraints you want in the syntax fmincon desires. The problem is that feasible sets consisting of disjoint, unconnected regions simply cannot be handled for theoretical reasons. fmincon has no way of jumping back and forth between the disjoint regions x(1)=0 and 0.4<=x(1)<=1 as it performs its search. Once fmincon lands in one of these regions, it will likely just stay bounded there and search for a local minimum.

Incidentally, you don't need three additional variables x(3), x(4), x(5) to do what you were attempting. You only need x(3) and the single additional constraint

 x(1)*(x(3)-x(1))=0

which is satisfied only when x(1) is equal either to 0 or to 0.4<=x(3)<=1.