fmincon for multiple variables

Hi,

I want to perform a constrained optimization,that makes fmincon a good option.

However, it returns the error '' fmincon requires all values returned by the functions to be of data type double''.

Y is the objective that I want to optimize,with three variables a1,a2,a3.

function Y = optfun(z)

a1 = z(1);

a2 = z(2);

a3 = z(3);

y1 = abs(vpa("0.129")*a1*sym(1i) + vpa("0.01")*a2*sym(1i) + vpa("8.57e-4")*a3*sym(1i) - vpa("2410.0i")/(vpa("50.0") + vpa("3.7e-5")*a3 + vpa("0.00111")*a2 + vpa("0.0333")*a1) - vpa("2410.0i")/(vpa("50.0") + vpa("2.96e-4")*a3 + vpa("0.00444")*a2 + vpa("0.0667")*a1) - vpa("2410.0i")/(vpa("50.0") + vpa("0.001")*a3 + vpa("0.01")*a2 + vpa("0.1")*a1) + vpa("25.0") + vpa("96.4i"))/abs(vpa("0.129")*a1*sym(1i) + vpa("0.01")*a2*sym(1i) + vpa("8.57e-4")*a3*sym(1i) + vpa("2410.0i")/(vpa("50.0") + vpa("3.7e-5")*a3 + vpa("0.00111")*a2 + vpa("0.0333")*a1) + vpa("2410.0i")/(vpa("50.0") + vpa("2.96e-4")*a3 + vpa("0.00444")*a2 + vpa("0.0667")*a1) + vpa("2410.0i")/(vpa("50.0") + vpa("0.001")*a3 + vpa("0.01")*a2 + vpa("0.1")*a1) + vpa("125.0") + vpa("96.4i"));

y2 = abs(vpa("0.173")*a1*sym(1i) + vpa("0.0135")*a2*sym(1i) + vpa("0.00115")*a3*sym(1i) - vpa("3250.0i")/(vpa("50.0") + vpa("3.7e-5")*a3 + vpa("0.00111")*a2 + vpa("0.0333")*a1) - vpa("3250.0i")/(vpa("50.0") + vpa("2.96e-4")*a3 + vpa("0.00444")*a2 + vpa("0.0667")*a1) - vpa("3250.0i")/(vpa("50.0") + vpa("0.001")*a3 + vpa("0.01")*a2 + vpa("0.1")*a1) + vpa("25.0") + vpa("130.0i"))/abs(vpa("0.173")*a1*sym(1i) + vpa("0.0135")*a2*sym(1i) + vpa("0.00115")*a3*sym(1i) + vpa("3250.0i")/(vpa("50.0") + vpa("3.7e-5")*a3 + vpa("0.00111")*a2 + vpa("0.0333")*a1) + vpa("3250.0i")/(vpa("50.0") + vpa("2.96e-4")*a3 + vpa("0.00444")*a2 + vpa("0.0667")*a1) + vpa("3250.0i")/(vpa("50.0") + vpa("0.001")*a3 + vpa("0.01")*a2 + vpa("0.1")*a1) + vpa("125.0") + vpa("130.0i"));

y3 = abs(vpa("0.197")*a1*sym(1i) + vpa("0.0153")*a2*sym(1i) + vpa("0.00131")*a3*sym(1i) - vpa("3690.0i")/(vpa("50.0") + vpa("3.7e-5")*a3 + vpa("0.00111")*a2 + vpa("0.0333")*a1) - vpa("3690.0i")/(vpa("50.0") + vpa("2.96e-4")*a3 + vpa("0.00444")*a2 + vpa("0.0667")*a1) - vpa("3690.0i")/(vpa("50.0") + vpa("0.001")*a3 + vpa("0.01")*a2 + vpa("0.1")*a1) + vpa("25.0") + vpa("148.0i"))/abs(vpa("0.197")*a1*sym(1i) + vpa("0.0153")*a2*sym(1i) + vpa("0.00131")*a3*sym(1i) + vpa("3690.0i")/(vpa("50.0") + vpa("3.7e-5")*a3 + vpa("0.00111")*a2 + vpa("0.0333")*a1) + vpa("3690.0i")/(vpa("50.0") + vpa("2.96e-4")*a3 + vpa("0.00444")*a2 + vpa("0.0667")*a1) + vpa("3690.0i")/(vpa("50.0") + vpa("0.001")*a3 + vpa("0.01")*a2 + vpa("0.1")*a1) + vpa("125.0") + vpa("148.0i"));

Y = (1-y1)+(1-y3)+(y2);

end

By the way,I want to know is there possibilities to add constration to my objective value?

For example, 0<[y1,y2,y3]<1.

Any help will be appreciate,thx!

[Stephen23]

Note: the context is explained in the comments here:

https://www.mathworks.com/matlabcentral/answers/2011672-optimization-of-multiple-variables

Answer 1

[Matt J]

Y = double( (1-y1)+(1-y3)+(y2) );

But I don't see why you are using vpa and other sym operations. I don't think it really increases precision for you since all your coefficients have substantially less than 32 decimal places.

[Matt J]

By the way,I want to know is there possibilities to add constration to my objective value? For example, 0<[y1,y2,y3]<1.

Of course! That's why it's called fmincon. For the constraints you describe, you would need to use the nonlcon argument.

[Guan Hao]

@Matt J The code works now!!!Thanks a lot :)

I'm using the optimization tool box like this:

But all the options I saw is to add constration on the variables a1 a2 a3.

I want to add them on y1 y2 y3.

[Matt J]

I want to add them on y1 y2 y3.

Yes, you will need a nonlinear constraint for that.

I have not used the optimization app, but presumably the "Nonlinear" button will prompt you to provide a constraint function handle as described here:

https://www.mathworks.com/help/optim/ug/fmincon.html#busog7r-nonlcon

[Guan Hao]

I think maybe I can create a constraint function like this:

function [c,ceq] = constraintFcn(z)
% Example:
% Constrain a solution to the region
% x^2 + y^2 <= 5
% x^2 + y^2 >= 2
% y = x^3

% Edit the lines below with your calculation
% Note, if no inequality constraints, specify c = []
% Note, if no equality constraints, specify ceq = []
y1 = z(1);
y2 = z(2);
y3 = z(3);
c(1) = y1-1;
c(2) = y3-1;
ceq=[];
end

[Matt J]

Yes, although this only implements 2 of the 6 bounds you said you wanted.

[Guan Hao]

@Matt J Hi,sorry to bother you again.

After I try my constraint function,the value didn't meet my requirement.(y1,y2,y3 returns the values of 0.2 0.2 0.2,respectively)

I'm trying to add constraint to my output function,not variables.

Is it correct to write the code like this?

function [c,ceq] = constraintFcn(z)

% My Constraint
% 0.9 <= y1 <= 1
% 0 <= y2 <= 0.1
% 0.9 <= y3 <= 1

y1 = z(1);
y2 = z(2);
y3 = z(3);
c(1) = -y1+0.9;
c(2) = y2-0.1;
c(3) = -y3+0.9;
c(4) = y1-1;
c(5) = y2;
c(6) = y3-1;
ceq=[];

end

[Torsten]

All these constraints can be set in the arrays "lb" and "ub" which are direct inputs to "fmincon":

lb = [0.9 0 0.9];
ub = [1 0.1 1]

And c(5) = y2 is incorrect because it sets y2 <= 0, not y2 >= 0.

[Guan Hao]

@Torsten

Thanks for your reply.

Isn't that setting bounds for the variables (a1 a2 a3)?

I want to set constrait to y1 y2 y3,which are functions of a1 a2 a3.

[Torsten]

Yes, you are right.

But c(5) should be -y2, not y2.

[Guan Hao]

@Torsten It returns the value of three variables [a1,a2,a3]=[0.9005,0.0126,0.9452],so I'm able to calculate my functions y1,y2,y3.

However, [y1,y2,y3]=[0.2,0.2,0.2].Obviously,they didin't follow my constraint for y1,y2,y3

[Torsten]

But you set the constraints on the parameters to be fitted, not on y1, y2 and y3 which you claim are functions of a1, a2 and a3. The vector z in "constraintFcn" is the vector (a1,a2,a3) from which you first have to compute y1,y2 and y3 (whatever they are). That's why I said that its the same as setting

lb = [0.9 0 0.9];
ub = [1 0.1 1]

And as I can see from your output: the vector x returned from "fmincon" satisfies lb <= x <= ub.

[Guan Hao]

@Torsten Okay...I see. But how to set the constraints to y1 y2 y3 ? In fact,I didn't need any constraints on a1 a2 a3,but there are boundary on y1 y2 y3

[Torsten]

Compute y1, y2 and y3 in "constraintFcn" from the z-vector and set the constraints as you did:

y = fun_z_to_y(z);
y1 = y(1);
y2 = y(2);
y3 = y(3);
c(1) = -y1+0.9;
c(2) = y2-0.1;
c(3) = -y3+0.9;
c(4) = y1-1;
c(5) = -y2;
c(6) = y3-1;

[Guan Hao]

@Torsten Thanks for your patience.

Did you mean something like this? It returns a error message that '' fmincon requires all values returned by the functions to be of data type double''.

Sorry,I must misunderstand your meaning...

function [c,ceq] = constraintFcn_new(z)

% My Constraint
% 0.9 <= y1 <= 1
% 0 <= y2 <= 0.1
% 0.9 <= y3 <= 1

a1 = z(1);

a2 = z(2);

a3 = z(3);

y = [abs(- vpa("0.122")*a1*sym(1i) - vpa("0.00895")*a2*sym(1i) - vpa("7.32e-4")*a3*sym(1i) + vpa("1530.0i")/(vpa("50.0") + vpa("0.001")*a3 + vpa("0.01")*a2 + vpa("0.1")*a1) + vpa("1530.0i")/(vpa("50.0") + vpa("8.0e-6")*a3 + vpa("4.0e-4")*a2 + vpa("0.02")*a1) + vpa("1530.0i")/(vpa("50.0") + vpa("6.4e-5")*a3 + vpa("0.0016")*a2 + vpa("0.04")*a1) + vpa("1530.0i")/(vpa("50.0") + vpa("5.12e-4")*a3 + vpa("0.0064")*a2 + vpa("0.08")*a1) + vpa("1530.0i")/(vpa("50.0") + vpa("2.16e-4")*a3 + vpa("0.0036")*a2 + vpa("0.06")*a1) - vpa("25.0") - vpa("102.0i"))/abs(vpa("0.122")*a1*sym(1i) + vpa("0.00895")*a2*sym(1i) + vpa("7.32e-4")*a3*sym(1i) + vpa("1530.0i")/(vpa("50.0") + vpa("0.001")*a3 + vpa("0.01")*a2 + vpa("0.1")*a1) + vpa("1530.0i")/(vpa("50.0") + vpa("8.0e-6")*a3 + vpa("4.0e-4")*a2 + vpa("0.02")*a1) + vpa("1530.0i")/(vpa("50.0") + vpa("6.4e-5")*a3 + vpa("0.0016")*a2 + vpa("0.04")*a1) + vpa("1530.0i")/(vpa("50.0") + vpa("5.12e-4")*a3 + vpa("0.0064")*a2 + vpa("0.08")*a1) + vpa("1530.0i")/(vpa("50.0") + vpa("2.16e-4")*a3 + vpa("0.0036")*a2 + vpa("0.06")*a1) + vpa("125.0") + vpa("102.0i"));

abs(- vpa("0.176")*a1*sym(1i) - vpa("0.0129")*a2*sym(1i) - vpa("0.00106")*a3*sym(1i) + vpa("2200.0i")/(vpa("50.0") + vpa("0.001")*a3 + vpa("0.01")*a2 + vpa("0.1")*a1) + vpa("2200.0i")/(vpa("50.0") + vpa("8.0e-6")*a3 + vpa("4.0e-4")*a2 + vpa("0.02")*a1) + vpa("2200.0i")/(vpa("50.0") + vpa("6.4e-5")*a3 + vpa("0.0016")*a2 + vpa("0.04")*a1) + vpa("2200.0i")/(vpa("50.0") + vpa("5.12e-4")*a3 + vpa("0.0064")*a2 + vpa("0.08")*a1) + vpa("2200.0i")/(vpa("50.0") + vpa("2.16e-4")*a3 + vpa("0.0036")*a2 + vpa("0.06")*a1) - vpa("25.0") - vpa("147.0i"))/abs(vpa("0.176")*a1*sym(1i) + vpa("0.0129")*a2*sym(1i) + vpa("0.00106")*a3*sym(1i) + vpa("2200.0i")/(vpa("50.0") + vpa("0.001")*a3 + vpa("0.01")*a2 + vpa("0.1")*a1) + vpa("2200.0i")/(vpa("50.0") + vpa("8.0e-6")*a3 + vpa("4.0e-4")*a2 + vpa("0.02")*a1) + vpa("2200.0i")/(vpa("50.0") + vpa("6.4e-5")*a3 + vpa("0.0016")*a2 + vpa("0.04")*a1) + vpa("2200.0i")/(vpa("50.0") + vpa("5.12e-4")*a3 + vpa("0.0064")*a2 + vpa("0.08")*a1) + vpa("2200.0i")/(vpa("50.0") + vpa("2.16e-4")*a3 + vpa("0.0036")*a2 + vpa("0.06")*a1) + vpa("125.0") + vpa("147.0i"));

abs(- vpa("0.223")*a1*sym(1i) - vpa("0.0163")*a2*sym(1i) - vpa("0.00134")*a3*sym(1i) + vpa("2790.0i")/(vpa("50.0") + vpa("0.001")*a3 + vpa("0.01")*a2 + vpa("0.1")*a1) + vpa("2790.0i")/(vpa("50.0") + vpa("8.0e-6")*a3 + vpa("4.0e-4")*a2 + vpa("0.02")*a1) + vpa("2790.0i")/(vpa("50.0") + vpa("6.4e-5")*a3 + vpa("0.0016")*a2 + vpa("0.04")*a1) + vpa("2790.0i")/(vpa("50.0") + vpa("5.12e-4")*a3 + vpa("0.0064")*a2 + vpa("0.08")*a1) + vpa("2790.0i")/(vpa("50.0") + vpa("2.16e-4")*a3 + vpa("0.0036")*a2 + vpa("0.06")*a1) - vpa("25.0") - vpa("186.0i"))/abs(vpa("0.223")*a1*sym(1i) + vpa("0.0163")*a2*sym(1i) + vpa("0.00134")*a3*sym(1i) + vpa("2790.0i")/(vpa("50.0") + vpa("0.001")*a3 + vpa("0.01")*a2 + vpa("0.1")*a1) + vpa("2790.0i")/(vpa("50.0") + vpa("8.0e-6")*a3 + vpa("4.0e-4")*a2 + vpa("0.02")*a1) + vpa("2790.0i")/(vpa("50.0") + vpa("6.4e-5")*a3 + vpa("0.0016")*a2 + vpa("0.04")*a1) + vpa("2790.0i")/(vpa("50.0") + vpa("5.12e-4")*a3 + vpa("0.0064")*a2 + vpa("0.08")*a1) + vpa("2790.0i")/(vpa("50.0") + vpa("2.16e-4")*a3 + vpa("0.0036")*a2 + vpa("0.06")*a1) + vpa("125.0") + vpa("186.0i"));];

y1 = y(1);
y2 = y(2);
y3 = y(3);
c(1) = -y1+0.9;
c(2) = y2-0.1;
c(3) = -y3+0.9;
c(4) = y1-1;
c(5) = -y2;
c(6) = y3-1;
ceq=[];

end

[Torsten]

Insert the line

y = double(y)

before you assign

y1 = y(1);
y2 = y(2);
y3 = y(3);
...
    

I wonder why you need symbolic variables to define y. Using numerical variables will usually speed up the numerical solver "fmincon" remarkably.

[Steven Lord]

Rather than duplicating your objective function's code in your constraint function, I'd define your objective function to return multiple outputs (the second and later ones being the quantities you're trying to constraint) and call it with multiple outputs in your constraint function.

Note that as per the description of the fun input argument on the fmincon documentation page, if you're using the SpecifyObjectiveGradient and HessianFcn options you may need to push these "supplementary outputs" later in the output list.

Y = optfun([1; 2; 3])

Y = 1.8024

confun([1; 2; 3])

  Y is 1.80239
 y1 is 0.198814
 y2 is 0.198692
 y3 is 0.197487
check: 1.80239

function [Y, y1, y2, y3] = optfun(z)
a1 = z(1);
a2 = z(2);
a3 = z(3);
y1 = abs(vpa("0.129")*a1*sym(1i) + vpa("0.01")*a2*sym(1i) + vpa("8.57e-4")*a3*sym(1i) - vpa("2410.0i")/(vpa("50.0") + vpa("3.7e-5")*a3 + vpa("0.00111")*a2 + vpa("0.0333")*a1) - vpa("2410.0i")/(vpa("50.0") + vpa("2.96e-4")*a3 + vpa("0.00444")*a2 + vpa("0.0667")*a1) - vpa("2410.0i")/(vpa("50.0") + vpa("0.001")*a3 + vpa("0.01")*a2 + vpa("0.1")*a1) + vpa("25.0") + vpa("96.4i"))/abs(vpa("0.129")*a1*sym(1i) + vpa("0.01")*a2*sym(1i) + vpa("8.57e-4")*a3*sym(1i) + vpa("2410.0i")/(vpa("50.0") + vpa("3.7e-5")*a3 + vpa("0.00111")*a2 + vpa("0.0333")*a1) + vpa("2410.0i")/(vpa("50.0") + vpa("2.96e-4")*a3 + vpa("0.00444")*a2 + vpa("0.0667")*a1) + vpa("2410.0i")/(vpa("50.0") + vpa("0.001")*a3 + vpa("0.01")*a2 + vpa("0.1")*a1) + vpa("125.0") + vpa("96.4i"));
y2 = abs(vpa("0.173")*a1*sym(1i) + vpa("0.0135")*a2*sym(1i) + vpa("0.00115")*a3*sym(1i) - vpa("3250.0i")/(vpa("50.0") + vpa("3.7e-5")*a3 + vpa("0.00111")*a2 + vpa("0.0333")*a1) - vpa("3250.0i")/(vpa("50.0") + vpa("2.96e-4")*a3 + vpa("0.00444")*a2 + vpa("0.0667")*a1) - vpa("3250.0i")/(vpa("50.0") + vpa("0.001")*a3 + vpa("0.01")*a2 + vpa("0.1")*a1) + vpa("25.0") + vpa("130.0i"))/abs(vpa("0.173")*a1*sym(1i) + vpa("0.0135")*a2*sym(1i) + vpa("0.00115")*a3*sym(1i) + vpa("3250.0i")/(vpa("50.0") + vpa("3.7e-5")*a3 + vpa("0.00111")*a2 + vpa("0.0333")*a1) + vpa("3250.0i")/(vpa("50.0") + vpa("2.96e-4")*a3 + vpa("0.00444")*a2 + vpa("0.0667")*a1) + vpa("3250.0i")/(vpa("50.0") + vpa("0.001")*a3 + vpa("0.01")*a2 + vpa("0.1")*a1) + vpa("125.0") + vpa("130.0i"));
y3 = abs(vpa("0.197")*a1*sym(1i) + vpa("0.0153")*a2*sym(1i) + vpa("0.00131")*a3*sym(1i) - vpa("3690.0i")/(vpa("50.0") + vpa("3.7e-5")*a3 + vpa("0.00111")*a2 + vpa("0.0333")*a1) - vpa("3690.0i")/(vpa("50.0") + vpa("2.96e-4")*a3 + vpa("0.00444")*a2 + vpa("0.0667")*a1) - vpa("3690.0i")/(vpa("50.0") + vpa("0.001")*a3 + vpa("0.01")*a2 + vpa("0.1")*a1) + vpa("25.0") + vpa("148.0i"))/abs(vpa("0.197")*a1*sym(1i) + vpa("0.0153")*a2*sym(1i) + vpa("0.00131")*a3*sym(1i) + vpa("3690.0i")/(vpa("50.0") + vpa("3.7e-5")*a3 + vpa("0.00111")*a2 + vpa("0.0333")*a1) + vpa("3690.0i")/(vpa("50.0") + vpa("2.96e-4")*a3 + vpa("0.00444")*a2 + vpa("0.0667")*a1) + vpa("3690.0i")/(vpa("50.0") + vpa("0.001")*a3 + vpa("0.01")*a2 + vpa("0.1")*a1) + vpa("125.0") + vpa("148.0i"));

% Converting to double after the fact is simpler to type for this example
% than removing all your vpa and sym calls, but for your real code I'd
% perform the removal to avoid involving Symbolic Math Toolbox at all
y1 = double(y1);
y2 = double(y2);
y3 = double(y3);
Y = (1-y1)+(1-y3)+(y2);
end

function confun(z)
% Reuse rather than duplicating
[Y, y1, y2, y3] = optfun(z);
fprintf("  Y is %g\n y1 is %g\n y2 is %g\n y3 is %g\n", Y, y1, y2, y3)
fprintf("check: %g\n", (1-y1)+(1-y3)+y2)
end

[Matt J]

Or better still, follow the guidelines here,

https://www.mathworks.com/help/optim/ug/objective-and-nonlinear-constraints-in-the-same-function.html

for avoiding repeated computations in the constraint and objective function code.

[Guan Hao]

@Steven Lord@Matt J

This is all of my code,it removes all vpa and sym things.

However,when I run the optimization,it seems that the constraints are on a1,a2,a3 again.

function [Y, y1, y2, y3] = optfun_new(z)

Z01=50;
Z02=75;
v=3e+8;
section=100;
L=0.1;
a=L/section;
a0=Z01;
f1=1e+9;
f2=2e+9;
w1=2*pi*f1;
w2=2*pi*f2;

syms w  

beta=w/v;

a1 = z(1);
a2 = z(2);
a3 = z(3);

% ABCD Matrix

for m=1:1:section

    l=m*a;

    Z0=a0+(a1*l)+(a2*(l^2))+(a3*(l^3));

    B(m)=1i*Z0*sin(beta*a);

    C(m)=1i*(1/Z0)*sin(beta*a);

end

Bsum=sum(B);

Csum=sum(C);

ABCD=[1 Bsum ; Csum 1];

% Convert into S11

S11=((ABCD(1,2)-ABCD(2,1)*Z01*Z02)+((ABCD(1,1)*Z02)-(ABCD(2,2)*Z01)))/((ABCD(1,2)+ABCD(2,1)*Z01*Z02)+((ABCD(1,1)*Z02)+(ABCD(2,2)*Z01)));

gamma=abs(S11);

% Optimize objective

y1=subs(gamma,w,w1);

y2=subs(gamma,w,((w1+w2)/2));

y3=subs(gamma,w,w2);

% Converting to double after the fact is simpler to type for this example
% than removing all your vpa and sym calls, but for your real code I'd
% perform the removal to avoid involving Symbolic Math Toolbox at all
y1 = double(y1);
y2 = double(y2);
y3 = double(y3);
Y = (1-y1)+(1-y3)+(y2);
end

function [c,ceq] = constraintFcn_new(z)
% Reuse rather than duplicating
[Y, y1, y2, y3] = optfun(z);
fprintf("  Y is %g\n y1 is %g\n y2 is %g\n y3 is %g\n", Y, y1, y2, y3)
fprintf("check: %g\n", (1-y1)+(1-y3)+y2)
y1 = z(1);
y2 = z(2);
y3 = z(3);
c(1) = -y1+0.9;
c(2) = y2-0.1;
c(3) = -y3+0.9;
c(4) = y1-1;
c(5) = -y2;
c(6) = y3-1;
ceq=[];
end

[Steven Lord]

You're calling the optfun to generate the intermediate values of your objective function, then after the fprintf lines (which I added just to show you that the constraint function had correctly obtained the values from the objective function) you throw away those intermediate values and replace them with elements of the input with which fmincon called the constraint function. Don't do that. Delete these lines from your constraint function:

y1 = z(1);
y2 = z(2);
y3 = z(3);

[Guan Hao]

@Steven Lord Thank you sir !!! It's working now.