Solving n x m non-linear equation systems where n > m

Hi all,

I am trying to solve m-nonlinear equations with n-unknowns. keeping in mind that n > m. What is the best way to do that using matlab?

My idea was to use "GlobalSearch" or "fsolve" functions in optimization toolbox matlab 2013b. the following is true for the problem in hand:

1. m = 2 & n = 4
2. f(x) = [@(x) (x(1)*cos(x(2)) - a*x(3) + b*x(4)) ; @(x) (x(1)*sin(x(2)) -b*x(3) + a*x(4))]
3. x(1),x(2),x(3),x(4) are the 4 unknowns variables in the equation
4. a,b are vectors of values that change with time.
5. The answer of those two equations is [y1,y2] which are vectors of values that change with time depending on a,b values.

from above this means that for variables a,b,y1,y2 at index "i" the equations will look like the follwoing:

f1(x) = y1(i) = x(1)*cos(x(2)) - a(i)*x(3) + b(i)*x(4)

f2(x) = y2(i) = x(1)*sin(x(2)) - b(i)*x(3) + a(i)*x(4)

the follwoing constraints applys:

• 0<= x(1) <= ((30e3*sqrt(2))+((30e3*sqrt(2))/2))
• -2*pi <= x(2) <= 2*pi
• 0<= x(3) <= inf
• -inf <= x(4) <= inf

rng default % For reproducibility
lb = [0,-2*pi,0,-inf];
ub = [30e3*sqrt(2),2*pi,inf,-inf];
X0 = rand(4,1); 
% Lets assume i = 500 and get the values at that index of our vectors
y1 = 2.4480e+04
y2 = -5.2580
a = -0.1333
b = 0.2775
Test.a = a;%OptimIn.a(iha);
Test.b = b;%OptimIn.b(iha);
Test.y1 = y1;%OptimIn.y1(iha);
Test.y2 = y2;%OptimIn.y2(iha);
GlobalSearch_Test.gs = GlobalSearch;
% for writing the Objective function i am not sure how to do it:
% I have the follwoing in general:
% Y = f(X,a,b); where Y = [y1;y2] and X = [x(1);x(2);x(3);x(4)]
% and i want to minimize the residual so:
% ß = min x sum(Y- f(X)); the sum is between i = 1 to i = end
GlobalSearch_Test.MyObj = @(x) ((x(1))^2 + (x(2))^2);% I am not sure what should the objective function be i want to minimize the residuales of the solution of those two equations
GlobalSearch_Test.c = @(x)[x(1)*cos(x(2)) - Test.a*x(3) + Test.b*x(4) - Test.y1;
            x(1)*sin(x(2)) - Test.b*x(3) + Test.a*x(4) - Test.y2];
GlobalSearch_Test.nonlinfcn = @(x)deal(GlobalSearch_Test.c(x));
GlobalSearch_Test.problem = createOptimProblem('fmincon','x0',X0,...
            'objective',GlobalSearch_Test.MyObj,'lb',lb,'ub',ub,...
            'nonlcon',GlobalSearch_Test.nonlinfcn);
GlobalSearch_Test.problem.options.TolCon = 1.0000e-50;
GlobalSearch_Test.x = run(GlobalSearch_Test.gs,GlobalSearch_Test.problem);
%% or the other method i tried is: 
% but i dont consider this the best way as i am not allowed to but limits for my search for optimal solution.
OptimIn.f = @(x) [(x(1)*cos(x(2)) - OptimIn.a*x(3) + OptimIn.b*x(4) - OptimIn.y1); (x(1)*sin(x(2)) - OptimIn.b*x(3) + OptimIn.a*x(4) - OptimIn.y2)];
OptimIn.X0 = rand(4,1);%zeros(4,1);%ones(4,1);
%OptimIn.Options = optimoptions('fsolve','Display','iter-detailed');
OptimIn.Options = optimset('Display','iter-detailed','TolFun',1e-50,'TolX',1e-50);
[Optim.X,Optim.fval,Optim.exitflag,Optim.output,Optim.jacobian] = fsolve(OptimIn.f,OptimIn.X0,OptimIn.Options);
    

so can you tell me what is the best way to solve such problem?

best regards

Answer 1

[Bjorn Gustavsson]

If you look at your equations, you can rewrite them (as far as I understand) to:

This you can treat as a linear system in x(3) and x(4) that you can solve explicitly then since the right-hand-side are variable with some limits on x(1) and x(2) this will correspond to 2 surfaces for x(3) and x(4). Here it might be best to explicitly calculate the inverse of your left-hand-side matrix just to multiply it with the right-hand-side-vector in order to get the explicit equations for the surfaces. Once you've done that you can evaluate that for all your values of interest. Since you want to calculate the surfaces for a -mesh it is easier (IMO) to do the explicit calculations of the matrix inverse and multiplication with the RHS - that way you get 2 simple expressions for x3 and x4 that you can easily calculate and display with for example pcolor.

HTH

[Aubai]

Hallo Bjorn,

Thanks for your fst replay. unfortionatily i dont understand the complete answer can you provide a simple code so i can understand what you are trying to say?

so what i think i understood from your answer is the following:

first i have to say my equations look something like that:

Y1 = x5.*cos(x6) -a.*x3 + b.*x4
Y2 = x5.*sin(x6) -b.*x3 + b.*x4

changing the equations will give us the following:

Y1 = X1 -a.*x3 + b.*x4
Y2 = X2 -b.*x3 + a.*x4

where

X1 = x5.*cos(x6)
X2 = x5.*sin(x6)

i changed the names of prvious x1 and x2 to make it easier to work with afterword.

then after i changed the equations as suggested i can solve it for x1,x2,x3,x4 using a very simple linear calcuation in matlab as following:

A = [1, 0, -a, b; 0, 1, -b, a];
b = [Y1; Y2];
X = A/b;

when looping for length(Y1) then we will get a solution for all available values of Y1, Y2, a,b in X which will have X1, X2 ,x3, x4 inside.

but using this method i am unable to give limits for x3 for example in which x3 should be always positive.

I am getting negative answer s for x3 how can i limit that to only positive soluations?

after that i am not sure how to understand or implement your suggestion

Best regards

[Bjorn Gustavsson]

OK, I missed the y1 and y2 components. But if you rewrite your equations with the non-linear terms on the RHS you should still get something like this:

Since you have a system of equations with larger number of unknowns you will get a set of solutions. Without the positivity-constraint on you would get full surfaces for and for the entire region of and . If you solve this system of equations you get:

For all values of x_1 and x_2 in the ranges and you will get some values for and , out of which you want to discard the regions where < 0. So lets do this:

X3 = @(a,b,y1,y2,x1,x2) ( (b*(y2 - x1.*sin(x2)))/(a^2 - b^2) - ...
                          (a*(y1 - x1.*cos(x2)))/(a^2 - b^2) );
X4 = @(a,b,y1,y2,x1,x2) ( (a*(y2 - x1.*sin(x2)))/(a^2 - b^2) - ...
                          (b*(y1 - x1.*cos(x2)))/(a^2 - b^2) ); 
% These functions should work for scalar parameters a b y1 and y2 with matrix-
% inputs for x1 and x2.
[x1,x2] = meshgrid(linspace(0,100,2001),linspace(-2*pi,2*pi,721));
% and the range of x1 is not your full range but you get the idea...
subplot(1,2,1)
pcolor(x1,x2,X3(1,2^.5,exp(3),log(4),x1,x2).*(X3(1,2^.5,exp(3),log(4),x1,x2)>=0)),shading flat % Same random selection of parameters
hold on
contour(x1,x2,X3(1,2^.5,exp(3),log(4),x1,x2),[0 0],'k')  
subplot(1,2,2)
pcolor(x1,x2,X4(1,2^.5,exp(3),log(4),x1,x2).*(X3(1,2^.5,exp(3),log(4),x1,x2)>=0)),shading flat % Same random selection of parameters
hold on
contour(x1,x2,X3(1,2^.5,exp(3),log(4),x1,x2),[0 0],'k')

That should give you the solutions for your x1, x2, x3 and x4, and you should also be able to see what regions give you solutions that fit your constraints. How you're going to select which solutions to use further you have to decide (or give additional information about...)

HTH

[Bjorn Gustavsson]

Did this give you enough of a solution?

[Aubai]

Still working on it that is why still no yes or no