objective and constraint functions definition using fmincon

Dear All,

Please advise me on this problem.

I have a vector of five unknown parameters (P=[P1 P2 P3 P4 P5]) and a matrix of N by 6 variables. (X=[Xij], i=1,...,N, j=1,..,6). (N is a large number).

The five unknown parameters are related nonlinearly to matrix entries X(i,j). As seen below, the const is a function of P vector, X matrix, where K and W are known constant matrices. ceq and X has the same dimensions.)

function [c,ceq] = const(P,X,W,K)   
c = [];
ceq(:,1) = K(1)*P(2)*X(:,1)-K(2)*(W(:,1)*sin(X(:,1))+(W(:,3)-W(:,2))*cos(X(:,1)))*(P(1)- 1)-W(:,5);
ceq(:,2) = K(1)*P(3)*X(:,2)-K(2)*sin(X(:,2))*(P(1) - 1)*W(:,1)-W(:,6);
ceq(:,3) = K(1)*P(4)*X(:,3)+.5*K(2)*(sin(X(:,3))*W(:,1)-(W(:,2)-W(:,3))*cos(X(:,3)))- W(:,5);
ceq(:,4) = K(1)*P(5)*X(:,4)+.5*K(2)*sin(X(:,4))-W(:,6);
ceq(:,5) = K(1)*P(2)*X(:,5)+K(2)*P(1)*(sin(X(:,5))*W(:,1)-(W(:,2)-W(:,3))*cos(X(:,5)))-W(:,5);
ceq(:,6) = K(1)*P(3)*X(:,6)+K(2)*P(1)*sin(X(:,6))-W(:,6);
end

The objective is to calculate the vector P (just five parameters) by minimizing the objective function obj. However, the obj(P1, X) is a function of just one of the unknown parameters and all the Xi,j variables.

 
obj =(25*(((2*P(1) - 1)*(cos(X(1,1))*cos(X(1,4))*sin(X(1,3))+……-cos(X(N,1))*cos(X(N,2))*(2*P(1)- 1))^2)^(1/2))/N);

In solving the optimization problem by fmincon function, I am not sure how to define the obj and how to call fmincon.

In defining the obj, what is the correct format? obj = @(?) (...)

Since we are only interested in finding the optimal values for P, how should I call the fmincon function? fmincon(obj,P0, [],[],[],[],[],[],@(?) const(P, X, W, K)) where P0 is the initial guess for the P vector.

Thanks in advance,

[Afsoon Nejati Aghdam]

I am adding more details to better ask for help.

function cost = parameterfun(P,X)
cost = 2*P(1) - 2*P(1)*(sin(X(1,5))*(cos(X(1,1))* ........(all the variables of X(i,j) but just P(1))..........sin(X(9,5))))^2)^(1/2);
end

After initialization of the following known constants:

W = [ Nx6]; K = [1x2]; lb = zeros(1,5); x0 = zeros(1,5); 

When I am running the following piece of code:

connL = @(P,X) const(P,X,W,K);
obj = @(P) parameterfun(P,X);
[Popt, fval] = fmincon(obj,x0,[],[],[],[],lb,[],connL);  

I receive the following error:

Unrecognized function or variable 'X'. Error in @(P)parameterfun(P,X) , Error in fmincon (line 566) initVals.f = feval(funfcn{3},X,varargin{:}); Caused by: Failure in initial objective function evaluation. FMINCON cannot continue.

Answer 1

[Walter Roberson]

W = [ Nx6]; K = [1x2]; lb = zeros(1,5); x0 = zeros(1,5); 

those initialize variables W and K (and lb and x0)

connL = @(P,X) const(P,X,W,K);

That creates an anonymous function . MATLAB sees the @(P,X) and takes note of the variable names listed there, P and X. It then scans the body of the anonymous function, and whereever it sees those variables it replaces them with references to parameters to be passed in, much like

connL = @(varargin) const(varargin{1}, varargin{2}, W, K)

As it encounters variables that are not on the parameter list, then at the time it builds the function it looks up the variables in scope, and copies them in to the function it is building, giving you something sort of like

global anon678735155_internal
anon678735155_internal.workspace{1}.W = W;    %copy as of time function is built
anon678735155_internal.workspace{1}.K = K;    %copy as of time function is built
anon678735155_internal.function = '@(P,X) const(P,X,W,K)';
mark_readonly anon678735155_internal
connL = @anon678735155

function output = anon678735155(varargin)
   global -readonly anon678735155_internal
   output = const(varargin{1}, varargin{2}, anon678735155_internal.workspace{1}.W, anon678735155_internal.workspace{1}.K);
end

(Not actual code, but works like this, if MATLAB offered a way to mark variables as read-only)

Then when connL is called upon, the values for W and K that were stored away are retrieved in order to execute the function call. And notice that assignments to W or K or P or X made in the body of the code after defining the anonymous function, have no effect on the execution of the anonymous function, which only relies on the the parameters passed in and the copies of the variables that it stored away.

obj = @(P) parameterfun(P,X);

That creates an anonymous function . MATLAB sees the @(P) and takes note of the variable name listed there, P. It then scans the body of the anonymous function, and whereever it sees that variable it replaces it with a eference to a parameter to be passed in, much like

obj = @(varargin) parameterfun(varargin{1},X)

As it encounters variables that are not on the parameter list, then at the time it builds the function it looks up the variables in scope, and copies them in to the function it is building, giving you something sort of like

global anon757740131_internal
anon757740131_internal.workspace{1}.X = X;    %copy as of time function is built
anon757740131_internal.function = '@(P) parameterfun(P,X)';
mark_readonly anon757740131_internal
obj = @anon757740131

function output = anon757740131(varargin)
   global -readonly anon757740131__internal
   output = parameterfun(varargin{1},anon678735155_internal.workspace{1}.X);
end

But -- there is no variable X in scope at the time that obj is being built, so it is not possible to store away a copy of the current X for later use.

When this happens, that the variable to be copied in is not defined in scope, then the effect is like

global anon757740131_internal
%no variables copied in
anon757740131_internal.function = '@(P) parameterfun(P,X)';
mark_readonly anon757740131_internal
obj = @anon757740131

function output = anon757740131(varargin)
   global -readonly anon757740131__internal
   output = parameterfun(varargin{1},anon678735155_internal.workspace{1}.X);
end

and then when you invoke obj and it gets to the point in execution that it needs X, it looks for anon678735155_internal.workspace{1}.X and sees that it does not exist, and complains about X not existing.

Remember that assignments in the code to the variable X made after the anonymous function is defined will not have an effect on the execution of the anonymous function: anonymous functions only rely on parameters passed in and stored copies of variables. In particular, if the anonymous function obj is later executed in a context that defined X then that X is ignored for this purpose. For example,

X = 5;
obj(73)

then when MATLAB executes anon757740131 it will not go looking in the calling environment to notice the X that was just defined. If a referenced variable in the body of an anonymous function does not exist then it gets marked as undefined, and the execution of the anonymous function fails as soon as the value is needed.

Please do not treat my code outline above as being exact: there are some complications having to do with global variables and shared variables.

[Popt, fval] = fmincon(obj,x0,[],[],[],[],lb,[],connL);  

When the non-linear constraint function, connL, is invoked, it will be passed the vector of current values. For example the first time, it would be invoked as connL(x0) . But you defined

connL = @(P,X) const(P,X,W,K);

which expects two parameters, not one.

If your X is known to your code, such as based on data read in, then you should define it in your code and then be using

X = appropriate value
W = [ Nx6]; K = [1x2]; lb = zeros(1,5); P0 = zeros(1,5); 

connL = @(P) const(P,X,W,K);
obj = @(P) parameterfun(P,X);
[Popt, fval] = fmincon(obj, P0, [], [], [], [], lb, [], connL);  

This then invokes the behaviour that the known X, W, K are copied into the anonymous function definitions . It differs from what you had in that X is defined ahead of time, and that X is no longer a parameter to the nonlinear constraint function connL. Also notice that I renamed x0 to P0 to emphasize that it is P that is being varied, not X.

[Afsoon Nejati Aghdam]

Dear Walter Roberson,

First, thank you so much for taking your time toward my question!

I'm not sure if I'm on the right track. Therefore, if there is something wrong, I would appreciate it if you could correct me.

The matrix X is related to vector P (both P and X are unknown) through the nonlinear equality constraints. For an initial guess of P0, I solved the nonlinear equations to obtain the equivalent initial guess of X0. (I used the fsolve function and checked the accuracy of the solution as well)

Then I defined the following variables and ran the following code.

W = [Nx6]; K = [1x2]; lb = zeros(1,5); P0 = [1x5];

X = X0;

connl = @(P) const(P,X,W,K);

obj = @(P) parameterfun(P,X);

[Popt, fval] = fmincon(obj,P0,[],[],[],[],lb,[],connl);

I receievd theis error:

Iter Func-count Fval Feasibility Step Length Norm of First-order

step optimality

0 6 3.649550e+01 1.914e+02 1.000e+00 0.000e+00 3.042e+01

1 12 6.724225e+01 6.098e-01 1.000e+00 1.615e+00 3.087e+05

2 68 6.724225e+01 6.098e-01 2.116e-09 7.356e-13 3.087e+05

Converged to an infeasible point.

fmincon stopped because the size of the current step is less than the value of the step size tolerance but constraints are not satisfied to within the value of the constraint tolerance.

<stopping criteria details>

Is the definition of X0 correct? In the above code, does matrix X vary during each iteration or considered constant (X0)?

Thanks so much in advance.

[Walter Roberson]

X would be considered constant, calculated once when you did the X = X0; and never calculated again.

[Walter Roberson]

If you need X to be calculated for each trial P value, then you need something like

W = [Nx6]; K = [1x2]; lb = zeros(1,5); P0 = [1x5]; 
connl = @(P) const(P,W,K);
obj = @parameterfun;
[Popt, fval] = fmincon(obj,P0,[],[],[],[],lb,[],connl);

function cost = parameterfun(P)
  X = calculate_X(P);
  cost = 2*P(1) - 2*P(1)*(sin(X(1,5))*(cos(X(1,1))* etc; %........(all the variables of X(i,j) but just P(1))..........sin(X(9,5))))^2)^(1/2);
end

function [c,ceq] = const(P,W,K)
    X = calculate_X(P);
    c = [];
    ceq(:,1) = K(1)*P(2)*X(:,1)-K(2)*(W(:,1)*sin(X(:,1))+(W(:,3)-W(:,2))*cos(X(:,1)))*(P(1)- 1)-W(:,5);
    ceq(:,2) = K(1)*P(3)*X(:,2)-K(2)*sin(X(:,2))*(P(1) - 1)*W(:,1)-W(:,6);
    ceq(:,3) = K(1)*P(4)*X(:,3)+.5*K(2)*(sin(X(:,3))*W(:,1)-(W(:,2)-W(:,3))*cos(X(:,3)))- W(:,5);
    ceq(:,4) = K(1)*P(5)*X(:,4)+.5*K(2)*sin(X(:,4))-W(:,6);
    ceq(:,5) = K(1)*P(2)*X(:,5)+K(2)*P(1)*(sin(X(:,5))*W(:,1)-(W(:,2)-W(:,3))*cos(X(:,5)))-W(:,5);
    ceq(:,6) = K(1)*P(3)*X(:,6)+K(2)*P(1)*sin(X(:,6))-W(:,6);
end

where calculate_X contains appropriate code.

If calculate_X is "expensive" then you could consider using memoize() https://www.mathworks.com/help/matlab/ref/memoize.html

It is tempting to think that you should calculate X in P and "global" it and have const pull out the value from global. However, it turns out that the order of calculation of the objective function and the nonlinear constraints is not fixed. With you having nonlinear equality constraints, most of the time fmincon() will probably call the nonlinear constraints several times in a row to estimate the gradient of the ceq entries to try to find a P that satisfies them all, and then only when it has satisfied the nonlinear constraints will it call the objective function. On the other hand, for the first several calls to the objective function, during which it is estimating the initial gradient, it ignores the nonlinear constraints. So constraints is not necessarily called immediately after objective, and objective is not necessarily called immediately after constraints.

[Afsoon Nejati Aghdam]

Dear Walter,

Thanks for the reply.

For the correct (optimal) value of vector P, the matrix X will be correct. So, for P0, we have X0, which is not the correct X. I do not need the values of X; I am interested in P only. However, since they are nonlinearly related, and at each trial, the P changes ... I guess the X should change accordingly. Finally, for the optimal value of P (Popt), we will have the correct value for matrix X (Xcorrect).

[Walter Roberson]

I guess the X should change accordingly

Then use the code structure I showed with both the objective function and the constraint function calling calculate_X

[Afsoon Nejati Aghdam]

Thanks for the comment. Just to add, I think using the connL is trivial. Since we calculate Xi for each trial from the nonlinear constraints, then there is no point in calling the connL again. It should be okay to go on with obj (no connL is needed!) and calling calculate_X inside obj.