Issues in the constraint functions of the Optimization Problem solving

I want to add more than two constraints which are linear, nonlinear and with both equality and inequality signs. As per the MATLAB documentation the examples shown are mostly of same variables which are written as array Examples.

I am using the following shown constraint for the optimisation. The constraint function is named as "area", this same function is called for the 'nonlcon' (default function for constraints).

#1 The equations commented are also to be considered as Constraints for the Optimisation. Do MATLAB optimisation considers this kind of constraint equations?

#2 How to check whether all constraints are satisfied.

function [c,ceq]= area(x0)
M=5;
x0=[5,8];

M2=m(M,x0(1));
M3=m(M2,x0(2));
M4=m(M3,(x0(1)+x0(2))); 

beta1=b(M,(x0(1)));
beta2=b(M2,(x0(2)));
beta3=b(M3,((x0(1)+x0(2))));

% M*sin(beta1)==M2*sin(beta2) & M2*sin(beta2)==M3*sin(beta3)
% (M4/M) <= 0.38
% beta1,beta2,beta3 < 62 degrees   (This can be written in terms of Thetas)
% Some more constraints will be added 

c=[];
ceq = log(2.4/((2.8*M*(sin(beta1)^2))-0.4)) + log(2.4/((2.8*M2*(sin(beta2)^2))-0.4)) + log(P(x0)) 

end

Below code shows the use of constraint function 'area' in the optimisation algorithm as "nonlcon".

%Initial Condition
% x0=[5,5] ;

% Lower bounds
lb=[1,1];


% Upper Bounds
ub=[30,40];

%constraint
nonlcon=@area

opts = optimoptions('fmincon','PlotFcn',["optimplotx","optimplotfunccount","optimplotfvalconstr","optimplotfval"],'Display','iter')
opts1 = optimoptions(opts,'MaxIterations',50,'StepTolerance',1e-9,'ConstraintTolerance',1e-9) 
[x,fval,exitflag,output,lambda,grad,hessian]= fmincon(@f,x0,[],[],[],[],lb,ub,nonlcon,opts1)

The optimisation is stopped between saying the below shown statement and then after enabling Feasibility Mode also the error is shown which is also mentioned in the end.

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>

Consider enabling the interior point method feasibility mode.


%Initial Condition
x0=[5,5] ;

% Lower bounds
lb=[1,1];


% Upper Bounds
ub=[30,40];

%constraint
nonlcon=@area
 

opts = optimoptions('fmincon','PlotFcn',["optimplotx","optimplotfunccount","optimplotfvalconstr","optimplotfval"],'Display','iter', ...
    'Algorithm','interior-point','EnableFeasibilityMode',true)
opts1 = optimoptions(opts,'MaxIterations',50,'StepTolerance',1e-9,'ConstraintTolerance',1e-9); % Recommended
[x,fval,exitflag,output,lambda,grad,hessian]= fmincon(@f,x0,[],[],[],[],lb,ub,nonlcon,opts1);




Converged to an infeasible point.

fmincon stopped because it is unable to find a point locally that satisfies
the constraints within the value of the constraint tolerance.

<stopping criteria details>

Even after changing the Constraint Tolerance it is not helping. Tried in the range 1e-9 to 1e-3

Thank you!

@Cris LaPierre @Walter Roberson

[Torsten]

What are m and b in the lines

M2=m(M,x0(1));
M3=m(M2,x0(2));
M4=m(M3,(x0(1)+x0(2))); 
beta1=b(M,(x0(1)));
beta2=b(M2,(x0(2)));
beta3=b(M3,((x0(1)+x0(2))));

?

[Walter Roberson]

beta1,beta2,beta3 < 62 degrees

do you need to be using sind()?

[Vivek]

@Torsten

The m and b are the external functions used to calculate the 'Mach' and 'Beta' for the given 'x0'

[Vivek]

@Walter Roberson These angles will be used in objective function and other constriants also . If we convert it to sind() the same range will be used in Objective function i.e., [-1,1]. So I guess it will change the values in Objective function.

Answer 1

[Sulaymon Eshkabilov]

There should be elementwise operations in ceq formulation if M2, beta1, beta2 and P are vectors.

Answer 2

[Matt J]

function [c,ceq]= area(x0)
M=5;
% x0=[5,8];

M2=m(M,x0(1));
M3=m(M2,x0(2));
M4=m(M3,(x0(1)+x0(2))); 

beta1=b(M,(x0(1)))*180/pi;
beta2=b(M2,(x0(2)))*180/pi;
beta3=b(M3,((x0(1)+x0(2))))*180/pi;

c(1)=M4-0.38*M;
c(2:4)=[beta1,beta2,beta3]-62;

ceq(1) = log(2.4/((2.8*M*(sin(beta1)^2))-0.4)) + log(2.4/((2.8*M2*(sin(beta2)^2))-0.4)) + log(P(x0)); 
ceq(2) =  M*sin(beta1)-M2*sin(beta2) ;
ceq(3) =  M2*sin(beta2)-M3*sin(beta3) ;
 
end

[Vivek]

@Matt J Thank you for the help!!

I tried this out and it is converging to Infeasible point.

Is there any way to check how many constraints are satisfied and how many did not?

[Vivek]

Also the equations c(1) and c(2) are <= given conditions.

When I try to change c(1)=M4-0.38*M; c(2:4)=[beta1,beta2,beta3]-62; to

c(1) >= M4-0.38*M;

c(2:4) >= [beta1,beta2,beta3]-62;

It is showing unrecognised 'c'

Also, I tried with using

c(1)= M4-0.38*M <=0;

c(2:4)= [beta1,beta2,beta3]-62 <=0;

This is showing error as

Error using fmincon

FMINCON requires all values returned by functions to be of data type double.

When used with ga() it is running properly and it is running as Unconstrained optimisation.

[Torsten]

An inequality constraint

 f(x) <= 0 

is specified as

c = f

not as

c = f<=0

An equality constraint

g(x)=0 

is specified as

ceq = g

not as

ceq = g==0 

[Vivek]

I tried with the below code, with the changes you had made. I am getting exitflag "-2" (I guess this means constraints are not satisfied).

function [c,ceq]= area(x0)
M=5;
% x0=[5,8];
M2=m(M,x0(1));
M3=m(M2,x0(2));
M4=m(M3,(x0(1)+x0(2))); 
beta1=b(M,(x0(1)))*180/pi;
beta2=b(M2,(x0(2)))*180/pi;
beta3=b(M3,((x0(1)+x0(2))))*180/pi;

c(1)=M4-0.38*M;
c(2:4)=[beta1,beta2,beta3]-62;
ceq(1) = log(2.4/((2.8*M*(sin(beta1)^2))-0.4)) + log(2.4/((2.8*M2*(sin(beta2)^2))-0.4)) + log(P(x0)); 
ceq(2) =  M*sin(beta1)-M2*sin(beta2) ;
ceq(3) =  M2*sin(beta2)-M3*sin(beta3) ;
 
end

The iterations ended with the message saying

"Optimization stopped because the relative changes in all elements of x are

less than options.StepTolerance = 1.000000e-08, but the relative maximum constraint

violation, 9.864905e-01, exceeds options.ConstraintTolerance = 1.000000e-05."

I am getting infeasible points.

@Matt J

[Torsten]

If you convert beta1, beta2 and beta3 in your code to degrees, I'm sure you had had to use

ceq(1) = log(2.4/((2.8*M*(sind(beta1)^2))-0.4)) + log(2.4/((2.8*M2*(sind(beta2)^2))-0.4)) + log(P(x0)); 
ceq(2) =  M*sind(beta1)-M2*sind(beta2) ;
ceq(3) =  M2*sind(beta2)-M3*sind(beta3) ;

[Vivek]

I was getting some real values which was infeasible with sin(beta). After using sind(beta) its returning complex values.

[Torsten]

Maybe you could try to take exp of your first equality constraint to avoid taking log of possibly negative expressions.

[Vivek]

Using exponential function will change the range of all numbers and it cant be used in further calculations.

[Matt J]

@Vivek I think you need to consider the possibility that your constraints might truly be infeasible. First, however, you should plot the constraint functions over the bounded region lb<=x<=ub and see if they are simultaneously satisfied anywhere. Because you only have two unknowns, this should be tractable. If you discover a region where they are satisfied, you should choose your initial guess x0 to lie in/near that region.

[Torsten]

Using exponential function will change the range of all numbers and it cant be used in further calculations.

Taking the exponential of a constraint does not change anything in the range of the parameters.

It doesn't matter whether you set

log(x)-2 = 0

or

exp(log(x)-2) = exp(0)

as a constraint.