incorporate the for - end loop in the objective function, where y is related with temp(x) rather than x

[mengyuli]

Dear Prof.,

If my objective function is min(y), where y=sum(temp(x)); 

x is the variable; temp here is the function of x, and can be calculated by a loop:

for i=1:length(x);

temp(i,:)=c(x(i));

end

My question is how I can incorporate the for - end loop in the objective function, where y is related with temp(x) rather than x? 

Thank you !

Best regards,

Mengyu

[Johan Löfberg]

Huh? I don't see the problem

I you already computed c, i.e it is a vector

y = sum(c)

or

y = 0;for i = 1:length(c);y = y + c(i);end

or if you mean c is a function

y = 0;for i = 1:length(x);y = y + c(x(i));end

or if c is a function which is vectorized

y = sum(c(x))

And I assume you mean that the objective is y, and y is a scalar, and you want to minimize it. In MATLAB min(y) is a function min which takes the minimum of the vector y, something completely different.

If not, you will have to explain better what you try to do

y = sum(c(x))

[mengyuli]

Thank you very much for your reply. My question actually is: 

x = sdpvar(10,1);

Constraints=[ x(1) == 1:5];

for i = 2:10;

    Constraints = [Constraints, x(i) == 1:8];

end

Objective=max(b(x));% I wanna minimise the maximum b(x), where b is defines by:

b=sum(c); %where c is a matrix and is the x-th row of the matrix d;

for i=1:10;

c(i,:)=d(x,1);

end

%solve the problem

ops = sdpsettings( 'solver','Cplex');

sol = optimize(Constraints,Objective,ops);

I tried to run the code, and an error occurs as ''Function 'subsindex' is not defined for values of class 'sdpvar'.''

I am guessing the error is because in the loop defining c,

for i=1:10;

c(i,:)=d(x,1);

end

So if I run the code directly, the Matlab will not recognize this code as  x here is the variable defined in x = sdpvar(10,1);

Do I need to define c as a constraint here?

Thank you very much.

[mengyuli]

I also tried to use another constraint to indicate c as:

x = sdpvar(10,1);

Constraints=[ x(1) == 1:5];

c=sdpvar(10,8760); 

Constraints=[ Constraints, implies(vx(1), c(1,:) == (squeeze(d(1,x(1),1,:)))']; %% d is a matrix with size of (1,49,10,8760), where c (size~(10,8760)) equals the last two dimentions of d;

for i = 2:10;

  Constraints = [Constraints, vx(i) == 1:num(1,i)];

  Constraints = [Constraints, implies(vx(i), c(i,:) == (squeeze(d(1,x(i),i,:)))'];

end

Objective=max(sum(c));

ops = sdpsettings( 'solver','Cplex');

sol = optimize(Constraints,Objective,ops);

However, an error occurs as ''Function 'subsindex' is not defined for values of class 'sdpvar'.''

 

If I don't use implies but use:  Constraints=[ Constraints, t(1,:) == (squeeze(temp0loadyear0(1,vx(1),1,:)))'];

An error occurs as ''Index exceeds matrix dimensions. Error in implies (line 50) Y = varargin{2};''.

[Johan Löfberg]

Code just doesn't make sense. x(1)  and x(i) etc are scalars. The data 1:5 and 1:8 are vectors of length 5 and 8. Hence you are saying the scalar variable x(i) should be a vector of some length. Absolutely impossible of course.

[mengyuli]

I wanna say the range of x(1) is a positive integer within the interval (1:5), i.e., x(1) can be 1, 2, 3, 4 or 5; 

[Johan Löfberg]

x = intvar(10,1);

Model = ismember(x(1),1:5)

[mengyuli]

This error occurs "Function 'subsindex' is not defined for values of class 'sdpvar'." when I wanna put x as an index in the matrix

t = sdpvar(1,49,200,8760);

Constraints = [Constraints, aux == temp0loadyear0];

for i=1:length(x);

    Constraints = [Constraints,t(i,:) == (squeeze(aux(1,x(i),i,:)))'];

end  

[Johan Löfberg]

you will have to supply reproducible code (as simple as possible still showing the errors)

[Johan Löfberg]

You cannot do that.

>> y = 1;

>> sdpvar x

>> y(x)

Function 'subsindex' is not defined for values of class 'sdpvar'.

You will have to come up with a logical model using binary variables and implies operator etc

[mengyuli]

Thank you very much for your time. Attached is my simplified code.

[Johan Löfberg]

well, it is solved already, you cannot index this way

[mengyuli]

Sorry, I didn't get it.  Is this :

d = binvar(10,1);

for i=1:length(vx);

    Constraints = [Constraints,implies(d(i),t(i,:)==(squeeze(aux(1,vx(i),i,:))))'];

end 

But I still have the same error. 

Thank you very much. 

[Johan Löfberg]

Instead of y = data(x(i)) you will have to do implies(x(i)==1, y == data(1), implies(x(i)==2, y == data(2)) etc.

[mengyuli]

Oh, ya! Thank you very much for your reply！ I have revised my code accordingly and it can be run now. However, the optimisation is extremely slow. I have attached in the attachment. Would you mind to see if the code can be optimised or not? I have tried Matlab global minimisation which seems faster. Is this normal? Thank you very much.

[mengyuli]

The size of temp0loadyear0.mat is 8MB, so I attached it in the Google drive link : 

https://drive.google.com/open?id=1-VOb8GVFuWtfaym_46zFQ7b-GkR6qcbv

[Johan Löfberg]

As I said, no need for any big model to reproduce this. You cannot index constant data using decision variables. You have to derive the logic model

[Johan Löfberg]

sorry, missed that you had fixed it

[Johan Löfberg]

A better model to do y = data(:,x) is

n = size(data,2);
d = binvar(n,1)
Model = [y == data*d, x == (1:n)*d, sum(d)==1]

[mengyuli]

Dear Prof., 

Thank you very much for your help. I have revised my constraints according to your suggestions and they work now. I am now trying to adjust my objective function as a Quadratic Programming, i.e.,

y1=demandeg(1:168)+t1./5-(sum(demandeg(1:168)+t1./5))/168;%where y1 is in the size of (168,1)

y2=sum(y1.*y1); % y2 should in the size of (1,1)

I am not sure if this problem belongs to Sum-of-squares programming. My code as attached works very slow.

The size of temp0loadyear1.mat is 11MB, so I also attached it in the Google drive link : 

https://drive.google.com/open?id=1wPUMfvlLmRY3VsOOEvN5tgySo4VDQWlR

Thank you very much!

[Johan Löfberg]

If you don't know for sure what sums-of-squares programming is, you definitely don't have a sums-of-squares program...

You have a mixed-integer quadratic program, so simply define your quadratic objective and use optimize as always.

[mengyuli]

Thank you very much for your reply. I tried to code my problem as a  quadratic objective, however, it takes a pretty long time to read the Objective1 function. 

%solve the problem

ops = sdpsettings( 'solver','Cplex');

Objective1=(sum(norm(demandeg(1:168)+t1./5-(sum(demandeg(1:168)+t1./5))/168),2))/168;

sol1 = optimize(Constraints,Objective1,ops);

vx2=value(vx);

opt1=value(Objective1);

Each time when running to the line defining Objective1, it is so slow. If I forced to quit the running, then the code gets stuck in the function y = times(X,Y). 

 May I know if there is another way that the problem can run faster if I still wanna use 'Cplex'? Many thanks.

[Johan Löfberg]

FYI, that is not a quadratic objective. That is a 2-norm objective. Equivalent to minimize, but not the same.

The sum operator is completely redundant. The norm of a vector is a scalar.

This should be slow and in theory impossible to solve. You have a mixed-integer second order cone program with more than 30000 variables...

https://yalmip.github.io/slowmilp

[Johan Löfberg]

Are you saying it gets stuck on Objective1=(sum(norm(demandeg(1:168)+t1./5-(sum(demandeg(1:168)+t1./5))/168),2))/168? 

That sounds weird. That does not happen here and should not be complecated for YALMIP to evaluate. Are you using the latest version?. In your previous code where you defined y2 as y2=sum(y1.*y1), that would definitely stall and take forever as you create a dense quadratic expression involving 900 million symbolic terms

[Johan Löfberg]

with mosek you have a pretty good solution rather quickly though, so reducing the optimality requirements you should be able to solve this (mosek finds a <1% suboptimal solution in 5 minutes). Haven't tested cplex

[mengyuli]

I tried to use mosek, and it got stuck at this stage for over half an hour:

K>>Objective1=(sum(norm(demandeg(1:168)+t1./5-(sum(demandeg(1:168)+t1./5))/168),2))/168;

K>> sol1 = optimize(Constraints,Objective1,ops);

Problem

  Name                   :                 

  Objective sense        : min             

  Type                   : CONIC (conic optimization problem)

  Constraints            : 45215           

  Cones                  : 1               

  Scalar variables       : 39493           

  Matrix variables       : 0               

  Integer variables      : 5723            

Optimizer started.

Mixed integer optimizer started.

Threads used: 4

Presolve started.

Presolve terminated. Time = 45.90

Presolved problem: 29708 variables, 24468 constraints, 4404857 non-zeros

Presolved problem: 0 general integer, 5438 binary, 24270 continuous

Clique table size: 207

BRANCHES RELAXS   ACT_NDS  DEPTH    BEST_INT_OBJ         BEST_RELAX_OBJ       REL_GAP(%)  TIME  

0        1        0        0        NA                   4.9988993063e-01     NA          60.0  

it seems it has return one result but still running.

[Johan Löfberg]

You've changed something, because the model you posted here has an optimal objective around 83, which also is the what the lower bound is in the root. That is not consistent with your display where a solution with value 0.4999 has been found. On my machine, I have a solution with objective 82.64 after 230 seconds, and a lower bound of 82.3.

But if you read the link I posted, you should know that this can take inf hours to terminate. You will have to play around with various solver options and tolerances and see if you can get something

[mengyuli]

This is the main code here. 

evnumberwd0=200;

vx = intvar(evnumberwd0,1);

t=sdpvar(168,evnumberwd0);

c=cell(evnumberwd0,1);

for i = 1:200; 

    c{i} = binvar(num(1,i),1);

end

i=1;

Constraints = [t(:,1) == (temp0loadyear1{1,1}(169:336,:))*c{1}, vx(1) == (1:num(1,1))*c{1}, sum(c{1})==1];

for i=2:200

Constraints = [Constraints,t(:,i) == (temp0loadyear1{1,i}(169:336,:))*c{i}, vx(i) == (1:num(1,i))*c{i}, sum(c{i})==1];

end

 

t1=sum(t,2);

ops = sdpsettings( 'solver','mosek');

%solve the problem

[Johan Löfberg]

Mosek solves it easily on my machine. Takes half an hour though, but compared to the age of the universe which is a very optimistic bound on the required time in theory, you should be ecstatic

Problem

  Name                   :                 

  Objective sense        : min             

  Type                   : CONIC (conic optimization problem)

  Constraints            : 45215           

  Cones                  : 1               

  Scalar variables       : 39493           

  Matrix variables       : 0               

  Integer variables      : 5723            

Optimizer started.

Mixed integer optimizer started.

Threads used: 2

Presolve started.

Presolve terminated. Time = 23.00

Presolved problem: 29708 variables, 24468 constraints, 4404857 non-zeros

Presolved problem: 0 general integer, 5438 binary, 24270 continuous

Clique table size: 207

BRANCHES RELAXS   ACT_NDS  DEPTH    BEST_INT_OBJ         BEST_RELAX_OBJ       REL_GAP(%)  TIME  

0        1        0        0        NA                   5.0073815427e-01     NA          32.7  

0        1        0        0        4.9241138202e-01     4.9241138202e-01     0.00e+00    1542.1

An optimal solution satisfying the relative gap tolerance of 1.00e-02(%) has been located.

The relative gap is 0.00e+00(%).

An optimal solution satisfying the absolute gap tolerance of 0.00e+00 has been located.

The absolute gap is 0.00e+00.

Objective of best integer solution : 4.924113820195e-01      

Best objective bound               : 5.007381542721e-01      

Construct solution objective       : Not employed

Construct solution # roundings     : 0

User objective cut value           : 0

Number of cuts generated           : 0

Number of branches                 : 0

Number of relaxations solved       : 1

Number of interior point iterations: 22

Number of simplex iterations       : 0

Time spend presolving the root     : 23.00

Time spend in the heuristic        : 0.00

Time spend in the sub optimizers   : 0.00

  Time spend optimizing the root   : 8.75

Mixed integer optimizer terminated. Time: 1542.17

Optimizer terminated. Time: 1542.31 

Integer solution solution summary

  Problem status  : PRIMAL_FEASIBLE

  Solution status : INTEGER_OPTIMAL

  Primal.  obj: 4.9241138202e-01    nrm: 8e+01    Viol.  con: 4e-13    var: 7e-16    cones: 0e+00    itg: 5e-14  

Optimizer summary

  Optimizer                 -                        time: 1542.31 

    Interior-point          - iterations : 0         time: 0.00    

      Basis identification  -                        time: 0.00    

        Primal              - iterations : 0         time: 0.00    

        Dual                - iterations : 0         time: 0.00    

        Clean primal        - iterations : 0         time: 0.00    

        Clean dual          - iterations : 0         time: 0.00    

    Simplex                 -                        time: 0.00    

      Primal simplex        - iterations : 0         time: 0.00    

      Dual simplex          - iterations : 0         time: 0.00    

    Mixed integer           - relaxations: 1         time: 1542.17 

[mengyuli]

Thank you very much for your reply. In my case, the solver run for a whole day and didn't stop. May I know if you have specified any parameters? 

[Johan Löfberg]

No. Mosek 8 no tricks, WIN64

[mengyuli]

So wired. Mine is Mac OS version

[Johan Löfberg]

Not really. Any epsilon change anywhere can cause an integer problem from being solved in no time to being impossible. Tha'ts the problem with nonconvex optimization.

You will have to test various solver options and see if it improves, or lower the requested tolerance to terminate early

(assuming that it hasn't simply crashed, i.e., that you actually see the branching process moving along)

[mengyuli]

thank you very much for your reply. If my variable number now is 2000 instead of 200, and my objective function now is linear as (Objective=max(dem(((1+(week-1)*168):(week*168)),ii)+t1./1000.*vari(ii))-max(dem((1+(week-1)*168):(week*168),ii))), where week varies from 1:52 and ii varies from 1:7. How to code properly to make code run faster if I have 52*7 objective functions, they are indepentent anyway.

vx = intvar(2000,1);

vx1=zeros(4,52,7,2000);

t=sdpvar (168,2000);

c=cell(2000,1);

ops = sdpsettings( 'solver','Cplex');

for i=1:4;

    c{1} = binvar(num(i,1),1);

    Constraints = [t(:,1) == (temp0loadweek1{1,i}(169:336,:))*c{1}, vx(1) == (1:num(i,1))*c{1}, sum(c{1})==1];

    for ii = 2:2000;

        c{ii} = binvar(num(i,ii),1);

        Constraints = [Constraints,t(:,ii) == (temp0loadweek1{ii,i}(169:336,:))*c{ii}, vx(ii) == (1:num(i,ii))*c{ii}, sum(c{ii})==1];

    end

    t1=sum(t,2);

    for week=1:52;

        for ii=1:7;

            Objective=max(dem(((1+(week-1)*168):(week*168)),ii)+t1./1000.*vari(ii))-max(dem((1+(week-1)*168):(week*168),ii));

            sol = optimize(Constraints,Objective,ops);

            vx1(i,week,ii,:)=value(vx);

        end

    end

end

[mengyuli]

I also tried to use parfor as:

for week=1:52;

        parfor ii=1:7;

            Objective=max(dem(((1+(week-1)*168):(week*168)),ii)+t1./1000.*vari(ii))-max(dem((1+(week-1)*168):(week*168),ii));

            sol = optimize(Constraints,Objective,ops);

            vx1(i,week,ii,:)=value(vx);

        end

end

Then there is an error as:

In sdpvar/saveobj (line 3)

  In parallel.internal.pool.serialize (line 21)

  In distcomp.remoteparfor (line 77)

  In parallel_function>iMakeRemoteParfor (line 1060)

  In parallel_function (line 444) 

Matrix dimensions must agree. 

It works with for -end loop, but is slow. 

[mengyuli]

The final goal of my code is to save optimised vx for 52 week in the vx1 tensor. 

[Johan Löfberg]

since you work with cells, there is no hope to make this significantly faster, as cell-based operations cannot be vectorized, and thus creates slow for-loop. If you would have had a matrix based format, you would be able to get rid of for-loops and indexing and write stuff like vx == integers*c, sum(c,1)==1 etc.

[Johan Löfberg]

parfor is not supported

[mengyuli]

vx == integers*c, sum(c,1)==1  This is for defining the constraints, and this for - loop runs for 1 min which is acceptable. 

For doing the optimisation, this for - loop takes a very long time. Is it possible to use optimizer here?

intvar week ii

Objective = max(dem(((1+(week-1)*168):(week*168)),ii)+t1./1000.*vari(ii))-max(dem((1+(week-1)*168):(week*168),ii));

 parfor week=1:52;

         for ii=1:length(vari);

            optimal=optimizer(Constraints, Objective,ops, {week , ii}, vx);

             vx1(i,week,ii,:)=optimal(week , ii);

         end

 end

However, there is an error saying subindex again is not supported for intvar

[Johan Löfberg]

Just s before you cannot do data(variable) where data is doubles and variable is an sdpvar. You have to model that manually using logic modelling and binary variables etc just as you've been told for your previous construct