How to use "fmincon" function to find the minimum of a function with nonlinear constraints

Everyone good time of day!

Please use the function "fmincon" to minimize the ellipse area function of the form , containing 6 variables from the canonical equation of the second order curve.

Everything would be nice if A, B, C, D, E and F were independent and linearly limited, but this is not so. These coefficients are limited by several inequalities of the form . Each such inequality suggests that the point with the coordinates is inside a certain second-order curve (in this case, an ellipse). The number of such inequalities corresponds to the length of the array of points, and usually at least twelve of them.

In addition to these inequalities, there are also two invariants that exclude the possibility of obtaining any other second-order curve except the ellipse: и .

In "MatLab help," I saw examples of using the function "fmincon" together with linear constraints in the form of inequalities, but I would like to put NONLINEAR constraints.

I will be incredibly grateful for the help!

P.S. I apologize for my English, because I am from Russia and do not know it perfectly.

Answer 1

[Alan Weiss]

I think that you are on the right track using fmincon. See the function reference page for details. Also, there is a topic on how to write nonlinear constraints.

Alan Weiss

MATLAB mathematical toolbox documentation

[Roman zuev]

Can you help me a little more?

I looked at the page that you offered, for which I thank you very much, but I was never able to "knock out" the desired result from the algorithm. As a result, it seems, an ellipse is built, but, at the same time, the wrong one - it is located incorrectly, and, most importantly, does not satisfy the main inequalities. The image I appended clearly shows that all points lie outside the perimeter, although the inequalities CLEARLY indicate that they ONLY NEED to be inside the perimeter of the ellipse they are looking for.

In the course of the "experiments," I noticed that the initial approximation - x0 - affects a lot and even in some cases (about this later) I get exactly what was required, but it remains a mystery to me why in some cases the answer is an ellipse that does not surround all points (for some reason), or even hyperbole (which cannot exist due to restrictions-invariants).

According to the console output in the command window, I can also assume that "fmincon" is looking for a local, not a global minimum, but why? What if there's a better option? And why is the minimum generally accepted that does not pass through the "barrier" of restrictions?

If it's not difficult for you, could you watch my code and somehow "fix" it by making it universal? I mean that, the code should work well at any array (not very large) of points at the input, and do not require manual adjustment of x0 values. After all, the fact is that this code is part of a large program, and the function "fmincon" will have to process various arrays of points several times.

Attached is a code in which, with the help of comments, I have tried to give as much detail as possible about what is happening.

P.S. I still failed to "force" the function to accept point arrays - XXX and YYY - in order not to re-write arrays within the function itself. For me, such a function challenge is something new, and my "old" methods do not work. I would be very grateful to you for your help with this:)

clear all

XXX = [15.78897133 55.27436793 43.437776    14.04490112  -23.10116586 -21.01148297]; % array X
YYY = [23.1372262  15.63395076 -27.45347755 -62.74204099 -67.99253449 -25.12781926]; % array Y

% A = x(1), B = x(2), C = x(3), D = x(4), E = x(5), F = x(6) - we rename as follows

% Sk0 = abs(-pi*(F-((D*(B*E - C*D))/(B^2 - A*C) - (E*(A*E - B*D))/(B^2 - A*C)))/(sqrt(A*C - B^2))); % ellipse area formula (exactly correct)
fun = @(x)abs(-pi*(x(6)-((x(4)*(x(2)*x(5) - x(3)*x(4)))/(x(2)^2 - x(1)*x(3)) - (x(5)*(x(1)*x(5) - x(2)*x(4)))/(x(2)^2 - x(1)*x(3))))/(sqrt(x(1)*x(3) - x(2)^2))); % ellipse area formula in new symbols
nonlcon = @ellipseparabola;
x0 = [1 2 3 4 5 6]; % 100 101 102 103 104 105 - works correct
A = []; % No other constraints
b = [];
Aeq = [];
beq = [];
lb = [];
ub = [];
% opts = optimoptions(@fmincon,'Algorithm','sqp'); % - - - -
sol = fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon);%, opts) %% this did not change anything

figure
hold on
axis equal
axis([min(XXX)-5 max(XXX)+5 min(YYY)-5 max(YYY)+5]);
plot(XXX,YYY,'o');
% A*x^2 + 2*B*x*y + C*y^2 + 2*D*x + 2*E*y + F = 0 - canonical equation of the second order curve
ellips_handle_function = @(x,y) sol(1)*x^2 + 2*sol(2)*x*y + sol(3)*y^2 + 2*sol(4)*x + 2*sol(5)*y + sol(6); % canonical equation of the second order curve in new symbols
fimplicit(ellips_handle_function);

function [c,ceq] = ellipseparabola(x)
    XXX = [15.78897133 55.27436793 43.437776    14.04490112  -23.10116586 -21.01148297]; % - - -
    YYY = [23.1372262  15.63395076 -27.45347755 -62.74204099 -67.99253449 -25.12781926]; % how can these lines be avoided?

    for i = 1:1:length(XXX) % for each point in the pattern, we create a constraint
        c(i) = x(1)*(XXX(i))^2 + 2*x(2)*XXX(i)*YYY(i) + x(3)*(YYY(i))^2 + 2*x(4)*XXX(i) + 2*x(5)*YYY(i) + x(6);% <= 0; if when you substitute a point inside the perimeter
        % in the canonical equation of a second-order curve, a negative number is obtained, then the point belongs to the perimeter - this is the condition and we write
    end
    
%     invar1 = A*C-B^2>0 or -A*C+B^2<0 :
    c(i+1) = -x(1)*x(3)+x(2)^2; % <= 0 - how to do simple "<" instead of "<=" ? this may be critical. (I tried to do this: if -A*C+B^2<=0, then -A*C+B^2 -0.0001 <0)
    delta = det([x(1) x(2) x(4); x(2) x(3) x(5); x(4) x(5) x(6)]); % matrix "delta" (part of the second invariant), the determinant of which is further
    I = x(1)+x(3); % some value "I" (part of the second invariant)
%     invar2 = delta*I<0;
    c(i+2) = delta*I; % <= 0 - not absolutely correctly
    ceq = [];
end

[Roman zuev]

Generally speaking, I found a solution to the problem I described, but if there is any other solution (better than mine) - I will be very glad: D