Non-linear system of 9 equations 9 unknowns MATLAB - unknowns coupled in polynomial ratio

Question

Goal

I want to solve a system of 9 non-linear equation with 9 unknowns with solve Matlab.

All 9 unknowns are coupled as a ration of polynoms (see myfun lower)

Fsolve

x02=[5000,5000,5000,0.4,0.4,0.4,0.4,0.4,0.4];
ctrl2=[9894+1i*0.118,9894+1i*0.118,9894+1i*0.118,0.5,0.5,0.5,0.5,0.5,0.5];
f2 = @(x) myfun(x,Peff);

options2 = optimoptions('fsolve','Algorithm','trust-region-dogleg','Display','iter-detailed'...
    ,'MaxFunEvals', 100000, 'MaxIter', 100000,'TolX',1e-12,'TolFun',1e-12,...
    'Jacobian','on');
[x2,F2,exitflag2,output2] = fsolve(f2,x02,options2);

Function myfun, return the system of equation to feed into fsolve an corresponding Jacobian

    function [F,J] = myfun( x, p)

% System of equation

F(1) = -(x(1) - x(1)*x(8)*x(9))/(x(4)*x(5) + x(6)*x(7) + x(8)*x(9) + x(4)*x(8)*x(7) + x(6)*x(5)*x(9) - 1) - p(1,1);
F(2) = -(x(2)*x(4) + x(2)*x(6)*x(9))/(x(4)*x(5) + x(6)*x(7) + x(8)*x(9) + x(4)*x(8)*x(7) + x(6)*x(5)*x(9) - 1) - p(1,2);
F(3) = -(x(3)*x(6) + x(3)*x(4)*x(8))/(x(4)*x(5) + x(6)*x(7) + x(8)*x(9) + x(4)*x(8)*x(7) + x(6)*x(5)*x(9) - 1) - p(1,3);
F(4) = -(x(1)*x(5) + x(1)*x(8)*x(7))/(x(4)*x(5) + x(6)*x(7) + x(8)*x(9) + x(4)*x(8)*x(7) + x(6)*x(5)*x(9) - 1) - p(2,1);
F(5) = -(x(2) - x(2)*x(6)*x(7))/(x(4)*x(5) + x(6)*x(7) + x(8)*x(9) + x(4)*x(8)*x(7) + x(6)*x(5)*x(9) - 1) - p(2,2);
F(6) = -(x(3)*x(8) + x(3)*x(6)*x(5))/(x(4)*x(5) + x(6)*x(7) + x(8)*x(9) + x(4)*x(8)*x(7) + x(6)*x(5)*x(9) - 1) - p(2,3);
F(7) = -(x(1)*x(7) + x(1)*x(5)*x(9))/(x(4)*x(5) + x(6)*x(7) + x(8)*x(9) + x(4)*x(8)*x(7) + x(6)*x(5)*x(9) - 1) - p(3,1);
F(8) = -(x(2)*x(9) + x(2)*x(4)*x(7))/(x(4)*x(5) + x(6)*x(7) + x(8)*x(9) + x(4)*x(8)*x(7) + x(6)*x(5)*x(9) - 1) - p(3,2);
F(9) = -(x(3) - x(3)*x(4)*x(5))/(x(4)*x(5) + x(6)*x(7) + x(8)*x(9) + x(4)*x(8)*x(7) + x(6)*x(5)*x(9) - 1) - p(3,3);

%% Jacobian

I compute the Jacobian myself but will spare you the detail, as it is considerably long

end

Results

                                         Norm of      First-order   Trust-region
 Iteration  Func-count     f(x)          step         optimality    radius
     0          1     2.45042e+19                      1.39e+14               1
     1          2     2.45042e+19              1       1.39e+14               1
     2          3     2.45031e+19           0.25       4.77e+15            0.25
     3          4     2.45031e+19          0.625       4.77e+15           0.625
     4          5     2.45031e+19        0.15625       4.77e+15           0.156
     5          6     2.44992e+19      0.0390625        6.8e+16          0.0391
     6          7     2.44992e+19      0.0976562        6.8e+16          0.0977
     7          8     2.44992e+19      0.0244141        6.8e+16          0.0244
     8          9      2.4495e+19     0.00610352       2.03e+17          0.0061
     9         10      2.4495e+19      0.0152588       2.03e+17          0.0153
    10         11      2.4486e+19      0.0038147       7.67e+17         0.00381
    11         12      2.4486e+19     0.00953674       7.67e+17         0.00954
    12         13     2.44592e+19     0.00238419       4.62e+18         0.00238
    13         14     2.44592e+19     0.00596046       4.62e+18         0.00596
    14         15     2.40048e+19     0.00149012       5.62e+20         0.00149
    15         16     2.40048e+19     0.00372529       5.62e+20         0.00373
    16         17     2.40048e+19    0.000931323       5.62e+20        0.000931
    17         18     2.40048e+19    0.000232831       5.62e+20        0.000233
    18         19     2.36832e+19    5.82077e-05       1.52e+21        5.82e-05
    19         20     2.36832e+19    0.000145519       1.52e+21        0.000146
    20         21      2.3131e+19    3.63798e-05       4.24e+21        3.64e-05
    21         22      2.3131e+19    9.09495e-05       4.24e+21        9.09e-05
    22         23     2.21355e+19    2.27374e-05       1.26e+22        2.27e-05
    23         24     2.21355e+19    5.68434e-05       1.26e+22        5.68e-05
    24         25     2.01772e+19    1.42109e-05        4.2e+22        1.42e-05
    25         26     2.01772e+19    3.55271e-05        4.2e+22        3.55e-05
    26         27      1.5592e+19    8.88178e-06       1.76e+23        8.88e-06
    27         28      1.5592e+19    2.22045e-05       1.76e+23        2.22e-05
    28         29     1.17854e+18    5.55112e-06       7.24e+23        5.55e-06
    29         30     3.43734e+16    1.38778e-05        1.9e+23        1.39e-05
    30         31     1.23843e+15    3.46945e-05       4.04e+22        3.47e-05
    31         32     1.23843e+15    8.67362e-05       4.04e+22        8.67e-05
    32         33     7.49991e+13     2.1684e-05       4.25e+21        2.17e-05
    33         34     7.49991e+13    5.42101e-05       4.25e+21        5.42e-05
    34         35     3.19073e+13    1.35525e-05       3.46e+21        1.36e-05
    35         36     3.19073e+13    3.38813e-05       3.46e+21        3.39e-05
    36         37     3.19073e+13    8.47033e-06       3.46e+21        8.47e-06
    37         38     3.19073e+13    2.11758e-06       3.46e+21        2.12e-06
    38         39     3.19073e+13    5.29396e-07       3.46e+21        5.29e-07
    39         40     3.19073e+13    1.32349e-07       3.46e+21        1.32e-07
    40         41     3.19073e+13    3.30872e-08       3.46e+21        3.31e-08
    41         42     3.19073e+13    8.27181e-09       3.46e+21        8.27e-09
    42         43     3.19073e+13    2.06795e-09       3.46e+21        2.07e-09
    43         44     3.19073e+13    5.16988e-10       3.46e+21        5.17e-10
    44         45     3.16764e+13    1.29247e-10          3e+21        1.29e-10
    45         46     3.16764e+13    1.29247e-10          3e+21        1.29e-10
    46         47     3.16764e+13    3.23117e-11          3e+21        3.23e-11
    47         48     3.16764e+13    8.07794e-12          3e+21        8.08e-12
    48         49     3.16764e+13    2.01948e-12          3e+21        2.02e-12
    49         50     3.16764e+13    5.04871e-13          3e+21        5.05e-13
    50         51     3.16764e+13    1.26218e-13          3e+21        1.26e-13
    51         52     3.16764e+13    3.15544e-14          3e+21        3.16e-14
    52         53     3.16764e+13    7.88861e-15          3e+21        7.89e-15
    53         54     3.16764e+13    1.97215e-15          3e+21        1.97e-15
    54         55     3.16764e+13    4.93038e-16          3e+21        4.93e-16

fsolve stopped because the relative norm of the current step, 1.097484e-16, is less than
max(options.TolX^2,eps) = 2.220446e-16. However, the sum of squared function values,
r = 3.167644e+13, exceeds sqrt(options.TolFun) = 1.000000e-06.

Optimization Metric                                               Options
relative norm(step) =   1.10e-16                max(TolX^2,eps) =   2e-16 (selected)
r =   3.17e+13                                    sqrt(TolFun) =  1.0e-06 (selected)


exitflag2 =

    -2


output2 = 

       iterations: 54
        funcCount: 55
        algorithm: 'trust-region-dogleg'
    firstorderopt: 3.000805388686251e+21
          message: 'No solution found.…'

Problem

Initial guess

x02 =

   1.0e+03 *

  Columns 1 through 5

   5.000000000000000   5.000000000000000   5.000000000000000   0.000400000000000   0.000400000000000

  Columns 6 through 9

   0.000400000000000   0.000400000000000   0.000400000000000   0.000400000000000

Solution given by fsolve

x2 =

   1.0e+03 *

  Columns 1 through 2

  5.000098340971978 - 0.000000066639557i  5.000100855522207 + 0.000000027141142i

  Columns 3 through 4

  5.000100887684736 + 0.000000021333305i  0.000500827051867 + 0.000000033172152i

  Columns 5 through 6

  0.000498312570833 - 0.000000060511167i  0.000500859436647 + 0.000000027409553i

  Columns 7 through 8

  0.000500831092720 + 0.000000033506374i  0.000500831171443 + 0.000000033543065i

  Column 9

  0.000498337065684 - 0.000000066909835i

What the solution should be (control)

    ctrl2 =

   1.0e+03 *

  Columns 1 through 2

  9.894000000000000 + 0.355900000000000i  9.894000000000000 + 0.355900000000000i

  Columns 3 through 4

  9.894000000000000 + 0.355900000000000i  0.000499999000000 + 0.000000000000000i

  Columns 5 through 6

  0.000499999000000 + 0.000000000000000i  0.000499999000000 + 0.000000000000000i

  Columns 7 through 8

  0.000499999000000 + 0.000000000000000i  0.000499999000000 + 0.000000000000000i

  Column 9

  0.000499999000000 + 0.000000000000000i

Comments

As you can see solve crashes fairly quick without going very far from the initial guess.

I tried changing TolX TolFun Algorithm but still crashes (following advice from Matlab website under what to do when algorithm fails)

The algorithm crashes similarly for other algorithm in fsolve and lsqnonlin (levnberg, trust-region-dogleg, trust-region-reflective)

Question to you:

Can you help figure out what I need to do in order to build an algorithm that will converge to the control values?

Thank you

Answer 1

What you are experiencing is actually nothing wrong with the solver. More than likely, you simply have many local minimizers in the problem, and the optimization algorithm is getting stuck there. There is really not much you can do about that... Almost all optimization solvers that you get your hand on are based on the principle of gradient descent (more or less), so they will guide you to a minimizer, but not necessarily the global one. What some people do is to generate a large set of random initial conditions, and to run the optimization solver for each initial condition. This isn't guaranteed to mean that one of the solutions you get back is a global minimizer, but at least you should be able to do better than a single initial guess.

If you really NEED to find the global minimizer, then you need to use some type of "Groebner-basis" solution method. These types of solvers are based on "elimination theory" and will yield a univariate polynomiial containing all of the global minimizers. However, these solvers are usually symbolic, extremely memory intensive, and basically aren't guaranteed to finish in a finite amount of time. I would suspect (I don't really know) that these problems are NP hard. In my experience, a seemingly small problem like the one you show can take very long to solve (if at all).