Error using fmincon function

*Hi everyone,

This is my first time using fmincon so I am sorry if the error is "stupid". Here is my code:*

 r1= sym('r1');  r2= sym('r2'); r3= sym('r3'); r4= sym('r4'); r5= sym('r5');

 syms A; syms B; syms C; syms alpha; syms beta; syms gamma;
 syms temp1; syms x;

 A = (r1^2)-(r2^2)-(r3^2)-(r4^2)+(2*r2*r4*sin(x/r5));
 B = (2*r2*r3*sin(x/r5))-(2*r2*r4);
 C = -2*r2*r3*cos(x/r5);

 alpha= (1+((B/C)^2));
 beta= (-2*A*B)/(C^2);
 gamma= ((A/C)^2)-1;

 temp1 = (-beta + (((beta^2)-(4*alpha*gamma))^(1/2)))/(2*alpha);

 fun = matlabFunction(temp1);

 rmax = 1;
 x0 = [0.5 0.5 0.5 0.5 0.5 0.5]';

 Matrix = [1 0 0 0 0 0;-1 0 0 0 0 0;0 1 0 0 0 0;0 -1 0 0 0 0;0 0 1 0 0 0;0 0 -1 0 0 0;0 0 0 1   0 0;0 0 0 -1 0 0;0 0 0 0 1 0;0 0 0 0 -1 0;0 0 0 0 0 1;0 0 0 0 0 -1];
 b = [rmax 0 0 0 0 0 0 0 0 0 0 0]';

 var = fmincon(fun,x0,Matrix,b);

Anyway, I am getting this error:

Warning: The default trust-region-reflective algorithm does not solve problems with the constraints you have specified. FMINCON will use the active-set algorithm instead. For information on applicable algorithms, see Choosing the Algorithm in the documentation. > In fmincon at 501 In Linkage_Optimization at 26 Warning: Your current settings will run a different algorithm (interior-point) in a future release.

> In fmincon at 506

In Linkage_Optimization at 26

Error using

symengine>makeFhandle/@(r1,r2,r3,r4,r5,x)(sqrt(-(1.0./r2.^2.*1.0./r3.^2.*1.0./cos(x./r5).^2.*(r2.*r4.*2.0-r2.*r3.*sin(x./r5).*2.0).^2+4.0).*(1.0./r2.^2.*1.0./r3.^2.*1.0./cos(x./r5).^2.*(-r1.^2+r2.^2+r3.^2+r4.^2-r2.*r4.*sin(x./r5).*2.0).^2.*(1.0./4.0)-1.0)+1.0./r2.^4.*1.0./r3.^4.*1.0./cos(x./r5).^4.*(r2.*r4.*2.0-r2.*r3.*sin(x./r5).*2.0).^2.*(r1.^2.*-2.0+r2.^2.*2.0+r3.^2.*2.0+r4.^2.*2.0-r2.*r4.*sin(x./r5).*4.0).^2.*(1.0./1.6e1))+1.0./r2.^2.*1.0./r3.^2.*1.0./cos(x./r5).^2.*(r2.*r4.*2.0-r2.*r3.*sin(x./r5).*2.0).*(r1.^2.*-2.0+r2.^2.*2.0+r3.^2.*2.0+r4.^2.*2.0-r2.*r4.*sin(x./r5).*4.0).*(1.0./4.0))./(1.0./r2.^2.*1.0./r3.^2.*1.0./cos(x./r5).^2.*(r2.*r4.*2.0-r2.*r3.*sin(x./r5).*2.0).^2.*(1.0./2.0)+2.0) Not enough input arguments.

Error in fmincon (line 640) initVals.f = feval(funfcn{3},X,varargin{:});

Error in Linkage_Optimization (line 26) var = fmincon(fun,x0,Matrix,b);

Caused by: Failure in initial user-supplied objective function evaluation. FMINCON cannot continue.

[Walter Roberson]

When you say

x0 = [0.5 0.5 0.5 0.5 0.5 0.5]'

which element of x0 corresponds to x, which to r1, r2, r3, r4, r5 ?

[Jonah]

I think the first should be r1, the second r2, the third r3, the fourth r4, the fifth r5, the sixth y.

Answer 1

[Walter Roberson]

fun = matlabFunction(temp1, 'vars', [r1, r2, r3, r4, r5, x]);
var = fmincon(@(X) fun(X(1),X(2),X(3),X(4),X(5),X(6)),x0,Matrix,b);

Or alternately, leave your fmincon function call alone and use

fun = matlabFunction(temp1, 'vars', {[r1 r2 r3 r4 r5 x]});

The syntax for this last version is not self-explanatory; look carefully at the Examples area of http://www.mathworks.com/help/symbolic/matlabfunction.html

[Jonah]

I am confused, I tried your first alternative but I get this error:

Warning: The default trust-region-reflective algorithm does not solve problems with the constraints you have specified. FMINCON will use the active-set algorithm instead. For information on applicable algorithms, see Choosing the Algorithm in the documentation. > In fmincon at 501 In Linkage_Optimization at 31 Warning: Your current settings will run a different algorithm (interior-point) in a future release. > In fmincon at 506 In Linkage_Optimization at 31

Solver stopped prematurely.

fmincon stopped because it exceeded the function evaluation limit, options.MaxFunEvals = 600 (the default value).

[Jonah]

Maybe a bit more information would help:

I am trying to maximize the last variable in my var vector while keeping the other variables smaller than rmax. Does that mean that I need to throw in a minus sign somewhere so that I find the maximum value instead of the minimum?

[Walter Roberson]

fun = matlabFunction(temp1, 'vars', [r1, r2, r3, r4, r5, x]);
options.MaxunEvals = 3000;   %more more more
options.Algorithm = 'active-set';
var = fmincon(@(X) fun(X(1),X(2),X(3),X(4),X(5),X(6)),x0,Matrix,b, options);

[Jonah]

Ok, I now get this error (when I change options to 3000):

Error using fmincon (line 226)
FMINCON requires the following inputs to be of data type double: 'Aeq'.

Error in Linkage_Optimization (line 32)
var = fmincon(@(X) fun(X(1),X(2),X(3),X(4),X(5),X(6)),x0,Matrix,b, options);

[Walter Roberson]

var = fmincon(@(X) fun(X(1),X(2),X(3),X(4),X(5),X(6)), x0, Matrix, b, [], [], options);

[Jonah]

So I tried this:

fun = matlabFunction(-temp1, 'vars', [r1, r2, r3, r4, r5, x]);

options.MaxFunEvals = 2799; 
Aeq = [];
beq = [];
lb = [];
ub = [];
nonlcon = [];

var = fmincon(@(X) fun(X(1),X(2),X(3),X(4),X(5),X(6)),x0,Matrix,b,Aeq,beq,lb,ub,nonlcon,options);

But I still get the same error:

    Solver stopped prematurely.

    fmincon stopped because it exceeded the function evaluation limit,
    options.MaxFunEvals = 2799 (the selected value).

[Jonah]

I feel like there is a deeper problem. I still need help. Any ideas? Here is my full current version of the code:

r1= sym('r1');  r2= sym('r2'); r3= sym('r3'); r4= sym('r4'); r5= sym('r5');

syms A; syms B; syms C; syms alpha; syms beta; syms gamma;
syms temp1; syms temp2; syms y1; syms y2; syms x;

A = (r1^2)-(r2^2)-(r3^2)-(r4^2)+(2*r2*r4*sin(x/r5));
B = (2*r2*r3*sin(x/r5))-(2*r2*r4);
C = -2*r2*r3*cos(x/r5);

alpha= (1+((B/C)^2));
beta= (-2*A*B)/(C^2);
gamma= ((A/C)^2)-1;

temp1 = (-beta + (((beta^2)-(4*alpha*gamma))^(1/2)))/(2*alpha);

rmax = 1;
x0 = [0.5 0.5 0.5 0.5 0.5 0.5]';

Matrix = [1 0 0 0 0 0;-1 0 0 0 0 0;0 1 0 0 0 0;0 -1 0 0 0 0;0 0 1 0 0 0;0 0 -1 0 0 0;0 0 0 1 0 0;0 0 0 -1 0 0;0 0 0 0 1 0;0 0 0 0 -1 0;0 0 0 0 0 1;0 0 0 0 0 -1];

b = [rmax 0 0 0 0 0 0 0 0 0 0 0]';

fun = matlabFunction(temp1, 'vars', [r1, r2, r3, r4, r5, x]);

options.MaxFunEvals = 2799; 
Aeq = [];
beq = [];
lb = [];
ub = [];
nonlcon = [];

var = fmincon(@(X) fun(X(1),X(2),X(3),X(4),X(5),X(6)),x0,Matrix,b,Aeq,beq,lb,ub,nonlcon,options);

[Walter Roberson]

Style note that would not affect anything: instead of the way you are using sym and syms now, you could use

syms r1 r2 r3 r4 r5 x

You do not need to mark A or B or alpha or beta or gamma or temp1 as syms: it will do that automatically because you are assigning symbolic expressions to those variables.

You are not currently using temp2 or y1 or y2

[Walter Roberson]

Your A, b bounds appear to be simple ones. I suggest that instead of using the A, b form, that you use lb and ub

Matrix = [];
b = [];
Aeq = [];
Beq = [];
lb = [0;    0; 0; 0; 0; 0];
ub = [rmax; 0; 0; 0; 0; 0];

Notice that you have constrained r2 to r5 and x to all be exactly 0 in your Matrix, b setup. If -1*x <= 0 and +1*x <= 0 then the only resolution is x = 0 exactly. Which would certainly go a ways towards explaining your difficulties. If variables are going to have inequalities strictly related to constants and not other variables, the lb / ub form is easier to read.

[Jonah]

Thanks for all the help thus far but I am still having issues. I am trying to optimize the lengths of r1, r2, r3, r4, and r5 so that the output is as large as possible. Essentially y = f(r1,r2,r3,r4,r5,x) and I am trying to find the maximum of y given a specific range for each r value. In actuality I am trying to model a mechanical linkage with each r being the length of a link in the linkage. The input is x and the output is y. I am trying to optimize the size of the links so that the output is as large as possible. The kinematic equation of the linkages is contained within my code. At this point the optimization seems to run but nothing happens, it just runs forever without finding a solution. I don't know what is wrong. Here is my current code based on your suggestions:

syms r1 r2 r3 r4 r5 x

A = (r1^2)-(r2^2)-(r3^2)-(r4^2)+(2*r2*r4*sin(x/r5));
B = (2*r2*r3*sin(x/r5))-(2*r2*r4);
C = -2*r2*r3*cos(x/r5);

alpha= (1+((B/C)^2));
beta= (-2*A*B)/(C^2);
gamma= ((A/C)^2)-1;

temp1 = (-beta + (((beta^2)-(4*alpha*gamma))^(1/2)))/(2*alpha);

fun = matlabFunction(temp1, 'vars', [r1, r2, r3, r4, r5, x]);

x0 = [1 2.5 2.5 2 0.5 1]';

options.MaxFunEvals = 100; 
Matrix = [];
b = [];
Aeq = [];
beq = [];
lb = [0.1; 0.1; 0.1; 0.1; 0.1; 0.1];
ub = [2; 4; 4; 4; 2; 5];
nonlcon = [];

var = fmincon(@(X) fun(X(1),X(2),X(3),X(4),X(5),X(6)),x0,Matrix,b,Aeq,beq,lb,ub,nonlcon,options);

[Walter Roberson]

One thing I would point out is that x ad r5 are only used in the pair x/r5 with neither of them used individually. You can therefor drop down a dimension in your search by using a single variable for the ratio.

With the constraints you have given just above for r5 and x, the ratio must be between (0.1/2) = 0.05 and (5/0.1) = 50.

However: You indicate "the input is x"; does that mean that in practice you will be given a particular x value, and need to optimize for that? If so then you really do need to drop down a dimension:

x = input('value for x?');
syms r1 r2 r3 r4 r5

A = (r1^2)-(r2^2)-(r3^2)-(r4^2)+(2*r2*r4*sin(x/r5));

B = (2*r2*r3*sin(x/r5))-(2*r2*r4);
C = -2*r2*r3*cos(x/r5);

alpha= (1+((B/C)^2));
beta= (-2*A*B)/(C^2);
gamma= ((A/C)^2)-1;

temp1 = (-beta + (((beta^2)-(4*alpha*gamma))^(1/2)))/(2*alpha);

fun = matlabFunction(temp1, 'vars', [r1, r2, r3, r4, r5]);

x0 = [1 2.5 2.5 2 0.5]';

options.MaxFunEvals = 100; 
Matrix = [];
b = [];
Aeq = [];
beq = [];
lb = [0.1; 0.1; 0.1; 0.1; 0.1];
ub = [2; 4; 4; 4; 2];

var = fmincon(@(R) fun(R(1),R(2),R(3),R(4),R(5)), x0, Matrix, b, Aeq, beq, lb, ub, nonlcon, options);

[Jonah]

Your reasoning seems to be sound. I tried your code but it does not output anything, it just gets stuck in fmincon.

[Walter Roberson]

Have you considered finding an approximate solution first? For example,

fragments = 16;
samples2 = linspace(0.1, 2, fragments+1);
samples4 = linspace(0.1, 4, fragments+1);
[R1, R2, R3, R4, R5] = ndgrid(samples2, samples4, samples4, samples4, samples2);

funvals = fun(R1, R2, R3, R4, R5);   %all of them at once

[minval, minidx] = min(funvals(:));
r1val = R1(minidx);
r2val = R2(minidx);
r3val = R3(minidx);
r4val = R4(minidx);
r5val = R5(minidx);

fprintf('minimum was %g at (%g, %g, %g, %g, %g)\n', minval, r1val, r2val, r3val, r4val, r5val);

After that you could take projections along various axes or do contour plots of iso slices, and so on, to figure out where the smaller values are clustered, and perhaps there are multiple such locations. Once the region of the lower solutions are found, you could zoom into that area by changing the ranges used in the ndgrid(), and then iterate again to get closer; and then take the ranges that you can extract from that and plop them into fmincon() to head for the minimum.

[Jonah]

First of all I am attaching an image that I just made to show the geometry of the linkage to give you a bit of physical significance (which always helps me conceptualize problems). How would I modify this code to find a maximum? I am trying to maximize the output. When I tried to simply change min(funvals(:)) to max, it gave me ridiculous outputs. What I am doing to check if the output is sane is using this simple code to plot it (I change the values in "fun" depending on the output of the code you provided):

syms r1 r2 r3 r4 r5 x

A = (r1^2)-(r2^2)-(r3^2)-(r4^2)+(2*r2*r4*sin(x/r5));
B = (2*r2*r3*sin(x/r5))-(2*r2*r4);
C = -2*r2*r3*cos(x/r5);

alpha= (1+((B/C)^2));
beta= (-2*A*B)/(C^2);
gamma= ((A/C)^2)-1;

temp1 = (-beta + (((beta^2)-(4*alpha*gamma))^(1/2)))/(2*alpha);

fun = matlabFunction(temp1, 'vars', [r1, r2, r3, r4, r5, x]);

x = 0:.1:10;

plot(x,fun(1,2.5,2.5,2,.5,x))
xlabel('x (input)');
ylabel('y (output)');

However I don't really even understand the physical significance of what I am plotting. Maybe it makes more sense to you? I assume that the best solution would be the maximum of output/input or in my case y/x for a given value of x. This is because I want my mechanical linkages to move a maximum amount for a given minimum input. If instance I move one link in the linkage 1cm (my input) and the output is 1.5cm. What do you think?

[Walter Roberson]

Remember, fmincon() always minimizes. If you have a function f() that you want to maximize, then minimize -1*f() .

[Jonah]

Right but what about in your "approximate solution"

[Walter Roberson]

The limit is indefinitely high.

r1 = smallest possible
r5 = 1
x = Pi/2
r3 = nearly largest possible
r4 same as r2
r2 = just slightly larger than r3, but not equal to r3

The closer r2 gets to r3, the higher the maximum

x/r5 = Pi/2 + n*2*Pi, n an integer

is a value at which singularities start to form. At that angle, the limit of f goes to (r1^2-(r2-r4)^2-r3^2) / (2*r2*(r3-r4))

you can drive this indefinitely high in more than one way; one of the ways involves maximizing r1 and ensuring that (r3-r4) is one side of 0, and the other involves minimizing r1 and ensuring that (r3-r4) is the other side of 0. And you can see from the numerator that you subtract off least if r4 approaches r2.

A different singularity forms along x/r5 = 3*Pi/2 + n*2*Pi, n integer. It can give you local maxima that do not depend upon the r3 vs r2 singularity. The form is very similar to the other:

-(r1^2-(r2+r4)^2-r3^2)/(2*r2*(r3+r4))

On the other hand, this can be driven indefinitely high by making r2 and r4 small, r1 small, and r3 high.

[Srikanth Gummadavalli]

@ walter-roberson

Hi

in the solution that you had provided

    fun = matlabFunction(temp1, 'vars', [r1, r2, r3, r4, r5, x]);
    var = fmincon(@(X) fun(X(1),X(2),X(3),X(4),X(5),X(6)),x0,Matrix,b);

what if we have r1 r2... up to r40. how to solve the above issue

Thankyou

[Walter Roberson]

r = sym('r', [1 40]);
fun = matlabFunction(temp1, 'vars', {r.'});
var = fmincon(fun, x0, Matrix, b)