FMINCON or other solvers???

[Ali_E]

hello professor, I want to optimize the following problem:

I've defined my problem in the following code:

clc;
clear;
close all;
%% Problem Definition
nVar=2;
VarSize=[1 nVar];
VarMin=10^(-3);
VarMax=10^2;
%% Weighted-Sum Approach
N=60;
w1=linspace(0,1,N);
w2=1-w1;
for i=1:N
s = [10^6 10^4 10^2];
S = diag(s);
deltat = sdpvar(3,3);
FWS=@(x) sum(w1(i)*Obj1(x)+w2(i)*Obj2(x)+deltat(i)'*S*deltat(i));
x0=unifrnd(VarMin,VarMax,VarSize);
A=[];
b=[];
LB=10^(-3);
UB=10^2;
options=optimoptions('fmincon','TolFun',10^(-12),'TolCon',10^(-10),'MaxFunEvals',4*10^3);
sol(i).Position=fmincon(FWS,x0,A,b,[],[],LB,UB,@NLC,options);
sol(i).Cost=MyCost(sol(i).Position);
end
Costs=[sol.Cost];
figure;
plot(Costs(1,:),Costs(2,:),'.');

I should say that deltat is a slack variable.

and I got error:
Supplied objective function must return a scalar value.

trying to debug I got that my objective function named FWS is:
Quadratic matrix variable 1x3

I should say that I defined Obj1 and Obj2 and MyCost and NLC as non linear constraints as follows:

function z = Obj1(x)
yref = 0;
p1 = -6.2039e-06;
p2 = 0.012729;
p3 = -0.040506;
p4 = -0.035217;
p5 = 0.28562;
p6 = -0.3329;
p7 = -0.19587;
p8 = 0.75971;
p9 = -0.54118;
p10 = -0.021845;
p11 = 0.098187;
mu = 50.5;
sigma = 29.011;
g = (x-mu)/sigma;
z = sum((yref-(p1*g.^10+p2*g.^9+p3*g.^8+p4*g.^7+p5*g.^6+p6*g.^5+p7*g.^4+p8*g.^3+p9*g.^2+p10*g+p11)).^2);
end

function z = Obj2(x)
yref = 0;
p1 = -0.0053913;
p2 = 0.011457;
p3 = 0.039991;
p4 = -0.1433;
p5 = 0.070996;
p6 = 0.31844;
p7 = -0.64456;
p8 = 0.33539;
p9 = 0.30767;
p10 = -0.37571;
p11 = 0.071788;
mu = 50.5;
sigma = 29.011;
g = (x-mu)/sigma;
z = sum((yref-(p1*g.^10+p2*g.^9+p3*g.^8+p4*g.^7+p5*g.^6+p6*g.^5+p7*g.^4+p8*g.^3+p9*g.^2+p10*g+p11)).^2);
end

function z = MyCost(x)
z1 = Obj1(x);
z2 = Obj2(x);
z = [z1;z2];
end

function [C, Ceq]=NLC(x)
deltat = sdpvar(3,3);
C=[((30767*((1000*x)/29011 - 50500/29011)^2)/100000 - (37571*x)/2901100 + (33539*((1000*x)/29011 - 50500/29011)^3)/100000 - (8057*((1000*x)/29011 - 50500/29011)^4)/12500 + (7961*((1000*x)/29011 - 50500/29011)^5)/25000 + (1278950236579183*((1000*x)/29011 - 50500/29011)^6)/18014398509481984 - (1433*((1000*x)/29011 - 50500/29011)^7)/10000 + (360206905396347*((1000*x)/29011 - 50500/29011)^8)/9007199254740992 + (6604510839140323*((1000*x)/29011 - 50500/29011)^9)/576460752303423488 - (3107872853893447*((1000*x)/29011 - 50500/29011)^10)/576460752303423488 + 37931111499167682390657/52261571515858183782400)^2 + ((787049070879267875*x)/1045231430317163675648 + (27059*((1000*x)/29011 - 50500/29011)^2)/50000 - (75971*((1000*x)/29011 - 50500/29011)^3)/100000 + (19587*((1000*x)/29011 - 50500/29011)^4)/100000 + (3329*((1000*x)/29011 - 50500/29011)^5)/10000 - (14281*((1000*x)/29011 - 50500/29011)^6)/50000 + (634413072308427*((1000*x)/29011 - 50500/29011)^7)/18014398509481984 + (2918764904100309*((1000*x)/29011 - 50500/29011)^8)/72057594037927936 - (3668884458035139*((1000*x)/29011 - 50500/29011)^9)/288230376151711744 + (1831068088942187*((1000*x)/29011 - 50500/29011)^10)/295147905179352825856 - 284748233055908747705/2090462860634327351296)^2-deltat
    -deltat];
Ceq=[];
end

how can I debug my program???

[Johan Löfberg]

fmincon expects the objective to be a scalar numerical value. You return a vector of sdpvar expressions?!?

>> FWS(x0)

Quadratic matrix variable 1x3 (full, real, 1 variable)

Coeffiecient range: 8.474 to 1000000

[Ali Esmaeilpour]

So which solver do you prefer for this nonlinear optimization?

[Ali Esmaeilpour]

How should I define slack variable without using sdpvar?

[Johan Löfberg]

I have absolutely no idea why you would think you require an sdpvar to introduce a slack when you are using low-level fmincon to solve a problem. 

[Johan Löfberg]

it looks like nasty nonlinear program, and you already setup a low-level code to use fmincon, so why not stay with fminon

[Ali Esmaeilpour]

It's because I dont Know How to define slack
Should I use syms?
or give an initial value to it and use there?

[Johan Löfberg]

you simply change your definition of x. first 2 elements are your original decision variables, and then the following are your slack decision variables. slack variables are not anything special, they are decision variables like any other decision variable in your model

[Ali Esmaeilpour]

So I Simply change those Obj as follows:
z=Obj1(x, delta)
.
.
.
z=funof(x) +delta'*s*delta

[Johan Löfberg]

no,  funof(x(1:n)) + x(n+1:end)'*s*x(n+1:end)

where n is the number of original variables

You anonymous function and fmincon assumes your objective is a function which is called with a vector of decision variables. In your case, that vector,i.e. the decision variables, is [originalvariables;slacks]

[Ali Esmaeilpour]

So we dont have delta and we have delta=x(n+1:end) as slack
But question there is if i got 2 original variables I'll have 4 slack variables because number of goals in my problem is 2. is this ok?

[Johan Löfberg]

I have no idea what you are doing so I cannot answer that.  The number of slacks/additional variables is simply the number of constraints you want to relax (with different slacks)

[Johan Löfberg]

from the picture delta has dimension w-1, hence the number of new variables is w-1 (although it is weird that you are using w functions in the objective, but then you only slack w-1 of them)

[Ali Esmaeilpour]

Another question is that I defined Obj1(x) and x is variable, but is it acceptable to write x(1:n)? Because x is a variable and nothing is pre-defined about x.

[Johan Löfberg]

Well you have to. You want Obj(oroginalvariables) + penalty(slacks). fmincon only sends a vector of decision variables, hence you have to have a way of partitioning the variable into the two 

[Ali Esmaeilpour]

Si I'll have:
z=Obj1(x)
.
.
.
n=2;
z=funof(x(1:n))+x(n+1:end)'*s*x(n+1:end)

It's weird because I haven't pre-defined x. and there is x(n+1:end) and What that "end" can mean?

[Johan Löfberg]

I am of course assuming still that you are using an anonymous function FWS=@(x) sum(w1(i)*Obj1(x(1:n)) + ...

Regarding "end", I can only advice you to study the basic of matlab

https://se.mathworks.com/company/newsletters/articles/matrix-indexing-in-matlab.html

[Ali Esmaeilpour]

I've already Sent you:

function z = Obj1(x)
yref = 0;
p1 = -6.2039e-06;
p2 = 0.012729;
p3 = -0.040506;
p4 = -0.035217;
p5 = 0.28562;
p6 = -0.3329;
p7 = -0.19587;
p8 = 0.75971;
p9 = -0.54118;
p10 = -0.021845;
p11 = 0.098187;
mu = 50.5;
sigma = 29.011;
g = (x-mu)/sigma;
z = sum((yref-(p1*g.^10+p2*g.^9+p3*g.^8+p4*g.^7+p5*g.^6+p6*g.^5+p7*g.^4+p8*g.^3+p9*g.^2+p10*g+p11)).^2);
end

and Obj2(x) is sth like above Function I've already Sent and these Obj form that FWS which you named anonymous

[Ali_E]

ok professor I changed all of those decision variables to what you suggested:

function z = Obj1(x)
yref = 0;
p1 = -6.2039e-06;
p2 = 0.012729;
p3 = -0.040506;
p4 = -0.035217;
p5 = 0.28562;
p6 = -0.3329;
p7 = -0.19587;
p8 = 0.75971;
p9 = -0.54118;
p10 = -0.021845;
p11 = 0.098187;
mu = 50.5;
sigma = 29.011;

s = [10^6 10^4 10^2];
S = diag(s);

n=3;
g = (x(1:n)-mu)/sigma;
z = sum((yref-(p1*g.^10+p2*g.^9+p3*g.^8+p4*g.^7+p5*g.^6+p6*g.^5+p7*g.^4+p8*g.^3+p9*g.^2+p10*g+p11)).^2)+x(n+1:end)'*S*x(n+1:end);
end

function z = Obj2(x)
yref = 0;
p1 = -0.0053913;
p2 = 0.011457;
p3 = 0.039991;
p4 = -0.1433;
p5 = 0.070996;
p6 = 0.31844;
p7 = -0.64456;
p8 = 0.33539;
p9 = 0.30767;
p10 = -0.37571;
p11 = 0.071788;
mu = 50.5;
sigma = 29.011;

s = [10^6 10^4 10^2];
S = diag(s);

n=3;
g = (x(1:n)-mu)/sigma;
z = sum((yref-(p1*g.^10+p2*g.^9+p3*g.^8+p4*g.^7+p5*g.^6+p6*g.^5+p7*g.^4+p8*g.^3+p9*g.^2+p10*g+p11)).^2)+x(n+1:end)'*S*x(n+1:end);
end

function [C, Ceq]=NLC(x)
n=3;
C=[((30767*((1000*x(1:n))/29011 - 50500/29011)^2)/100000 - (37571*x(1:n))/2901100 + (33539*((1000*x(1:n))/29011 - 50500/29011)^3)/100000 - (8057*((1000*x(1:n))/29011 - 50500/29011)^4)/12500 + (7961*((1000*x(1:n))/29011 - 50500/29011)^5)/25000 + (1278950236579183*((1000*x(1:n))/29011 - 50500/29011)^6)/18014398509481984 - (1433*((1000*x(1:n))/29011 - 50500/29011)^7)/10000 + (360206905396347*((1000*x(1:n))/29011 - 50500/29011)^8)/9007199254740992 + (6604510839140323*((1000*x(1:n))/29011 - 50500/29011)^9)/576460752303423488 - (3107872853893447*((1000*x(1:n))/29011 - 50500/29011)^10)/576460752303423488 + 37931111499167682390657/52261571515858183782400)^2 + ((787049070879267875*x(1:n))/1045231430317163675648 + (27059*((1000*x(1:n))/29011 - 50500/29011)^2)/50000 - (75971*((1000*x(1:n))/29011 - 50500/29011)^3)/100000 + (19587*((1000*x(1:n))/29011 - 50500/29011)^4)/100000 + (3329*((1000*x(1:n))/29011 - 50500/29011)^5)/10000 - (14281*((1000*x(1:n))/29011 - 50500/29011)^6)/50000 + (634413072308427*((1000*x(1:n))/29011 - 50500/29011)^7)/18014398509481984 + (2918764904100309*((1000*x(1:n))/29011 - 50500/29011)^8)/72057594037927936 - (3668884458035139*((1000*x(1:n))/29011 - 50500/29011)^9)/288230376151711744 + (1831068088942187*((1000*x(1:n))/29011 - 50500/29011)^10)/295147905179352825856 - 284748233055908747705/2090462860634327351296)^2-x(n+1:end)
    -x(n+1:end)];
Ceq=[];
end

function z = MyCost(x)
z1 = Obj1(x);
z2 = Obj2(x);
z = [z1;z2];
end

clc;
clear;
close all;
%% Problem Definition

nVar=4;

VarSize=[1 nVar];
VarMin=10^(-3);
VarMax=10^2;
%% Weighted-Sum Approach
N=60;
w1=linspace(0,1,N);
w2=1-w1;
for i=1:

N
FWS=@(x) w1(i)*Obj1(x)+w2(i)*Obj2(x);

x0=unifrnd(VarMin,VarMax,VarSize);
A=[];
b=[];
LB=10^(-3);
UB=10^2;
options=optimoptions('fmincon','TolFun',10^(-12),'TolCon',10^(-10),'MaxFunEvals',4*10^3);
sol(i).Position=fmincon(FWS,x0,A,b,[],[],LB,UB,@NLC,options);
sol(i).Cost=MyCost(sol(i).Position);
end
Costs=[sol.Cost];
figure;
plot(Costs(1,:),Costs(2,:),'.');

but FWS is a matrix now and it says it has to be a scalar!!!

[Johan Löfberg]

Of course it does as you've defined it as such...

[Johan Löfberg]

or maybe this is what you meant- Even worse

>> w1(i)*Obj1(x0)

ans =

     0     0     0

     0     0     0

     0     0     0

>> w2(i)*Obj2(x0)

ans =

   1.0e+09 *

    9.3101    0.0000    0.0000

    0.0000    0.0931    0.0000

    0.0000    0.0000    0.0009

This is turning into a support thread for basic MATLAB, which we cannot have here.