How to define integral function in objective function of YALMIP code

[Nur]

Hi,

I need to solve mixed integer optimization problem having objective function which consists several parts like quadratic function, linear function and most critical two integral functions. I am using MOSEK as external solver. The mathematical model is given in the attached PDF file.

I am in dire need to know whether it is possible to define integral function in objective function in direct or indirect  way …if so , Can you please tell me how ????

I am currently using code for unit commitment example given in YALMIP wiki examples.

http://users.isy.liu.se/johanl/yalmip/pmwiki.php?n=Examples.UnitCommitment

Thanks in advance,

Regards,

Nur

[Johan Löfberg]

Not possible (in neither Mosek nor YALMIP). You would have to come up with, e.g., a piecewise affine approximation of the integral.

[Johan Löfberg]

BTW, your notation is very strange, making it hard to see if this perhaps is trivial.

You write

alpha(Pw-P)

and then you write

alpha int (w-wi)f(w)dw

Is alpha a constant or a function? (looks like function first, but then outside integral, so constant). What are the decision variables and what are constants? Is the only decision variable in the integral w_si. If so, that whole expression simplifies to

-w_si*alpha*integral + alpha*integral

i.e, linear in w_si, simply compute the integral and then plug it into the optimization problem

[Nur]

Hi Johan,

Thanks for the reply. sorry for the ambiguity in the expression. I am trying to rephrase my optimization problem for clarity and posting it again (Unit Commitment reposted.pdf) along with the code (main file UC_main.m along with other function, xls files).

It would be great help if any suggestion is made to solve this sort of  optimization problem in any ways (directly or with approximation) using YALMIP platform.

Please suggest.  I am in dire need.

Thanks again.

Regards,

Nur

 

 

[Johan Löfberg]

Still not clear, as you have a subindex W on the function f. Does the function in the integrand depend on the decision variable

I.e, with y being the decision variable, do you have

f(y) = int_0^y (y-x)f(x)dx

or

f(y) = int_0^y (y-x)f(x,y)dx 

[Nur]

sorry Johan for making things difficult. In my case, the function f_w  only depends on integrand variable, f(x)dx  (not on the decision variable,y ) which is the first case of your given example.

My case,

f(y) = int_0^y (y-x)f(x)dx

[Johan Löfberg]

Actually it doesn't make any difference. Anyhow, unless there is something trivial about the integrals allowing you to compute them analytically and deriving a linear or quadratic expressions, you have to approximate them. Since you seem to want to stay in MILP/MIQP-workd, which is good, I would suggest approximating the integrals, as a function of decision variable, using a SOS2 PWA approximation

http://users.isy.liu.se/johanl/yalmip/pmwiki.php?n=Commands.SOS2

All you need is function values of the integrals, which easily are computed, e.g., (here with f = sin^3, decision variable allowed to take values between 0 and 2*pi, hence we need samples there)

Q = @(P,y) (y-P).*sin(P).^3;

values = (0:0.01:2*pi);

for i = 1:length(values)

    f1(i) = integral(@(P)Q(P,values(i)),0,values(i)); 

    f2(i) = integral(@(P)Q(P,values(i)),values(i),values(end)); 

end

plot(values,[f1' f2'])

[Nur]

Thanks a lot Johan for your quick support. I will now try to implement your suggestion in the code.

[Nur]

Hi Johan!!

I have implemented the integral function of the objective function (of unit commitment) using SOS2 PWA approximation according  to your suggestion. Now I have minimized objective function and found out the result also.

Is there a possibility that due to the SOS2 PWA approximation and the nature of objective function (which can be non-convex in nature , I am not sure), I am getting local optima not global optima?

How to find it out? 

Any suggestion/insight will be appreciated.

Thanks in advance.

Regards 

Nur

Note: I am now currently using GUROBI since it supports the SOS2 PWA approximation. I have attached the UC mathematical model and necessary codes again for your convenience.

[Johan Löfberg]

It is global up to the precision of the approximation you've selected (and optimality tolerances in the solver of course)

[Johan Löfberg]

Is this really correct

    Constraints = [Constraints, P_w(:,k) == lambda1(:,k)'*P_w_p', y1(k) == lambda1(:,k)'*f1(k,:)',lambda1(:,k)>=0,sum(lambda1(:,k))==1];
    Constraints = [Constraints, P_w(:,k) == lambda2(:,k)'*P_w_p', y2(k) == lambda2(:,k)'*f2(k,:)',lambda2(:,k)>=0,sum(lambda2(:,k))==1];

You have P_w(:,k) on the left-hand side in both constraints

The reason I am doubting the model in general is the curious fact that all lambda terms are 0,1, or 0.71 or 0.28 (when I used P_w_p=linspace(0,wr,10), similiar when using other levels of discretization)

[Nur]

Thanks Johan for looking into the code.

In the constraints what I tried to do is to implement it in the same way as:

http://users.isy.liu.se/johanl/yalmip/pmwiki.php?n=Commands.SOS2

In the objective function the 2nd (linear part),  3rd and 4th term(integral terms) has same decision variable  P_w. This is the reason I have

put same variable P_w in both of the constraints and in the linear part of objective function with an intention that while minimizing P_w it takes

into account all the three terms.

Objective = Objective + P(:,k)'*Cp*P(:,k) + Cq'*P(:,k)+Cr'*OO(:,k)+(S_w**P_w(:,k)*)+(y1(k)+y2(k));

After seeing your observation, I looked into the values of y1, y2 besides lambda1/lambda2. Though the values of lambda1 and lambda2 are same,the values of y1 and y2 are different.

I am confused about the implication of same values for both lambdas. Is it ok for my case (the 2nd (linear part),  3rd and 4th term (integral terms) has same decision variable  P_w. )?



Thanks

Nur

[Johan Löfberg]

Sorry, see now what you do, y1 and y2 both represent functions of P_w

[Nur]

Dear Johan,

I am back with another problem.  This is the continuation and extension of previous problem. Previously I was facing a problem with solving integral function (having decision variable in limits ) in the objective function. You came up with solution of specially ordered set type 2 (SOS2). It worked fine. Thanks for the solution. 

Now I am trying to integrate EV charging scheduling optimization in the same objective function using the difference between two integral functions (last term of objective function in the attached pdf). I tried to implement that in my objective function but It is showing no suitable solver. Currently I am using gurobi solver along with Yalmip. Is it possible to this in any other way??

Objective = Objective + P(:,k)'*Cp*P(:,k) + Cq'*P(:,k)+Cr'*OO(:,k)+(S_w'*P_w(:,k))+z(:,k)'*(y1(:,k)+y2(:,k))+((sum(y1(:,k),1))-(sum(y2(:,k),1)))*ones(1,N_clev)*((P_ev(:,k)));

It would be really helpful if you can give a possible solution of this problem.

Thanks in advance.

Nur

N.B.: I have attached PDF showing objective function and constraints and codes with the mail. I have successfully used the code for optimization upto fourth terms. Only problem is arising when I tried to integrate 5th term :(

[Johan Löfberg]

Looks to me as if you are trying to multiply two continuous decision variables. That is not possible with MILP/MIQP solvers as it is a nonconvex term. Not fun to deal with. The only approach I can think of of is triangulate the domain of interest, and then use a piecewise affine approximation in each simplex. No fun at all to setup.

[Johan Löfberg]

Well, a bit fun...I had to try and see. The code below is very naive (uniform gridding, and simple generation of MILP model. There is a whole field devoted to research on how to create these pwa models. I recommend various papers by Juan Pablo Vielma)

[X,Y] = meshgrid(-1:0.01:1,-1:0.01:1);
% Global minimum 0 at [0.5;-0.25]
mesh(X,Y,sin((X-.5).^2+(Y+.25).^2).^2)

% Create a triangulation
x = sdpvar(2,1);
f = @(x) sin((x(1)-.5)^2+(x(2)+.25)^2)^2;
% 2n simplicies will be generated
n = 5;
[f_pwa,model] = pwatriangulate2D([-1 <= x <= 1],f,x,n);

% And try to find the solution
solvesdp(model,f_pwa)

% Actual value in found point vs approximation value
 [f(double(x)) double(f_pwa)]

function [f_pwa,model] = pwatriangulate2D(D,f,x,n);

[~,L,U] = boundingbox(D,[],x);

X = linspace(L(1),U(1),n);
Y = linspace(L(2),U(2),n);

% Number of squares (two simplicies in each)
in_region = binvar((length(X)-1)*(length(Y)-1),1);
k = 1;

f_pwa = 0;
% We are in one box (save bounds also)
model = [D, sum(in_region) == 1];
for i = 1:length(X)-1
    for j = 1:length(Y)-1
        % Two simplicies in each box
        which_simplex = binvar(2,1);
        v1 = [X(i);Y(j)];
        v2 = [X(i+1);Y(j)];
        v3 = [X(i);Y(j+1)];
        v4 = [X(i+1);Y(j+1)];
        z1 = f(v1);
        z2 = f(v2);
        z3 = f(v3);
        z4 = f(v4);
        % Interpolation variables in two simplicies
        corners1 = sdpvar(3,1);
        corners2 = sdpvar(3,1);
        x1 = sdpvar(2,1);
        x2 = sdpvar(2,1);
        model = [model, implies(in_region(k), x == x1+x2), corners1 >=0, corners2 >= 0];
       
        % Interpolate function and variable in first simplex
        f_pwa = f_pwa + z1*corners1(1)+z2*corners1(2)+z3*corners1(3);
        model = [model, x1 == v1*corners1(1)+v2*corners1(2)+v3*corners1(3)];                
        model = [model, sum(corners1) == which_simplex(1)];        
       
        % Second simplex
        f_pwa = f_pwa + z2*corners2(1)+z3*corners2(2)+z4*corners2(3);
        model = [model, x2 == v2*corners2(1)+v3*corners2(2)+v4*corners2(3)];                
        model = [model, sum(corners2) == which_simplex(2)];
               
        % We are in at most one simplex, if we are in this region
        model = [model, sum(which_simplex) == in_region(k)];
        k = k+1;      
    end
end        

[Johan Löfberg]

A better approach if you use a MIQP solver could be to write x*y as .5*(x+y)^2-.5x^2-.5*y^2. The first quadratic term is convex and can thus be added to the objective function. The two nonconvex terms can be handled using sos2 approximations. In other words, instead of using one nasty 2D pwa which will involve O(m^2 ) binaries where m is a measure of the approximation level, you get two simple SOS2 models each involving O(m) binaries etc. If some term in the product is modeled using a SOS2 model already, you should of course reuse the intermediate variables involved in that model.

[Johan Löfberg]

like this

[X,Y] = meshgrid(-1:0.01:1,-1:0.01:1);

ftemp = cos(X+sin(X)).^2;
% f(x) + f(x)*y = f(x) + .5(f(x)+y)^2-.5*f(x)^2-.5*y^2
mesh(X,Y,ftemp + ftemp*Y);

sdpvar x y f f2 y2
n = 25;
xyi = linspace(-1,1,n)';
fi = cos(xyi+sin(xyi)).^2;
fi_squared = fi.^2;
lambda1 = sdpvar(length(xyi),1); % for f and f^2
lambda2 = sdpvar(length(xyi),1); % for y^2
F = [sos2(lambda1), sos2(lambda2)];
F = [F, x == lambda1'*xyi, f == lambda1'*fi, lambda1>=0,sum(lambda1)==1]
F = [F, f2 == lambda1'*fi_squared];
F = [F, y == lambda2'*xyi, y2 == lambda2'*xyi.^2, lambda2>=0,sum(lambda2)==1]

Objective = f + 0.5*(f+y)^2-0.5*f2-0.5*y2;
solvesdp([-1 <= [x;y] <= 1,F],Objective)

[Nur]

Thanks Johan for the quick reply ... I am looking into this.

[Nur]

Dear Johan,

I have implemented your suggestion according to your previous example. It is not converging to acceptable limit even after running for hours. Rather it is stuck at high gap of around 16 %.

I dont have any idea what may be responsible for this occurance and any possible way out of this problem.

I have attached PDF file explaing the problem with necessary equations (page-3).

Please help me with your valuable suggestions.

Regards,

Nur

[Johan Löfberg]

Well, we've gone from a very very very hard problem to a very very hard problem. To arrive at a just very hard problem, much more clever ideas are probably required.

What strikes me though is your extremely fine and uniform gridding in your pwa approximations. You should definitely try to reduce those (no reason to use lots of points in regions where the function is almost flat (hessian close to zero)

Tweaking solver settings can also be helpful. Sometimes small changes give huge impact (in both ways...)

[Johan Löfberg]

You should really take a look at the data you have.

These f2 functions, which you model with different sos2 pwa models, are exactly the same functions w.r.t the third argument

plot([f2(3,:,1)' f2(3,:,2)' f2(3,:,3)'])

and they simply scale linearly w.r.t to the first argument

plot([f2(3,:,1)'./f2(2,:,1)' f2(1,:,2)'./f2(2,:,2)'])

Hence, all your f2 functions can be represented using a single pwa model

the same holds for the f1 function