nonlinear equations system

morbidlizard · **Joined:** Sat Jun 04, 2011 3:58 pm **Posts:** 5

I have a system of non-linear equations which has to be solved.
This problem was previously solved by means of mathcad, but now it's considered unconvinient and I have to program an application with cute GUI etc. But, the problem is, I can't find any way to solve the system.
Mathcad claims it uses levenberq-marquardt method to solve my system, so, obviously, I tried to use Alglib's minlm* routines. But the algorithm stops because of "relative step is no more than EpsX"(terminationtype == 2), and the solution provided is not as good as I want it to be (and not as Mathcad's solution is).
So I'm looking for any advice on how to get a proper result (or which other algorithm to use).
Thanks in advance.

Sergey.Bochkanov · **Joined:** Fri May 07, 2010 7:06 am **Posts:** 927

Can you post your equation here along with starting point? And, if you can, your code which calls ALGLIB. We can try to find the reason for the failure to converge.

I am really interested in solving this problem because I want ALGLIB to be as good as possible.

morbidlizard · **Joined:** Sat Jun 04, 2011 3:58 pm **Posts:** 5

the original system is in .png file attached to this reply

Attachment:

File comment: system

system.png [ 22.22 KiB | Viewed 12213 times ]

unknown variables are L, H, A and R_x. other variables will be taken from the input but now I treat them as constants for testing purpose.
the code is as follows

Code:

using System;
using System.Text;

namespace Project1
{
    class Program
    {
        //init
        static double d = 3;
        static double a = 1000;
        static double F_0 = 0.035;
        static double F_H = 549;
        static double F_d = 9;
        static double c_yd = 1.3;
        static double c_xd = 0.5;
        static double F_sh = 5;
        static double c_ysh = 0.6;
        static double c_xsh = 0.9;
        static double M_g = 1500;
        static double M_z = 400;
        static double hamma_M = 74000;
        static double V = 2.2;
        //
        static double ro = 1030;
        static double v = 1.3e-6;
        static double g = 9.80665;
        static double d_a = d / a;
        static double F_f = F_H * a / d;
        static double G_g = M_g * g * (hamma_M - ro * g) / hamma_M;
        static double G_z = M_z * g * (hamma_M - ro * g) / hamma_M;
        static double R_yd = c_yd * ro * V * V * F_d / 2;
        static double R_xd = c_xd * ro * V * V * F_d / 2;
        static double Q = c_ysh * ro * V * V * F_sh / 2;
        static double R_xsh = c_xsh * ro * V * V * F_sh / 2;
        
        //
        static void f(double[] x, double[] fi, object obj)
        {
            double L = x[0];
            double H = x[1];
            double A = x[2];
            double R_x = x[3];
            fi[0] = R_x - (0.04 * (A + 110.0 * Math.Exp(-6.74 * F_0) * (Math.PI * H * L / (4.0 * F_f) - 0.02)) - 0.09) * (ro * V * V * F_H / 2.0);

            fi[1] = Math.PI * H * L / (4.0 * F_f) - 1.0 + Math.Exp(0.0008 * (2.0 * R_yd + Q + 2.0 * G_g + G_z) / ((R_x + R_xsh) * F_0));

            fi[2] = A + 110.0 * Math.Exp(-6.74 * F_0) * (Math.PI * H * L / (4.0 * F_f) - 0.02) - A * (1.0 - Math.Exp(-4.5 * (2.0 * R_yd + Q + 2.0 * G_g + G_z) / (R_x + R_xsh)));

            fi[3] = H / L - Math.Exp(1.0 - 2.0 * R_yd / (Q + 2.0 * G_g + G_z));
        }

        static void Main()
        {




            double[] x = { 50, 50, 6, 10000 };

            double[] fi = { 0, 0, 0, 0 };

            double epsg = 0.0000000000001;
            double epsf = 0;
            double epsx = 0;
            int maxits = 0;
            alglib.minlmstate state;
            alglib.minlmreport rep;

            alglib.minlmcreatev(4, x, 0.0001, out state);
            alglib.minlmsetcond(state, epsg, epsf, epsx, maxits);
            alglib.minlmoptimize(state, f, null, null);
            alglib.minlmresults(state, out x, out rep);

            System.Console.WriteLine("{0}", rep.terminationtype); // EXPECTED: 4
            
            f(x, fi, null);
            System.Console.WriteLine(" X : {0}, {1}, {2}, {3}", x[0], x[1], x[2], x[3]);
            System.Console.WriteLine(" f(X) : {0}, {1}, {2}, {3}", fi[0], fi[1], fi[2], fi[3]);
            System.Console.ReadLine();
        }
    }
}

it's almost completely copied from alglib example for minlv*, save for the function and dimension chanched
I sincerely hope, there is any hint on what to do
Anyway, thank you for your responsiveness

Sergey.Bochkanov · **Joined:** Fri May 07, 2010 7:06 am **Posts:** 927

I suggest you to scale your variables with alglib.minlmsetscale() function. Appropriate scale can be [10, 10, 10, 10000]. Unfortunately, I've been too busy yesterday to try it out, but I hope that I will find some time today.

Sergey.Bochkanov · **Joined:** Fri May 07, 2010 7:06 am **Posts:** 927

I've tested your code. The problem was caused by different scale of the variables, which: a) caused numerical differentiation to lose precision, b) caused incorrect work of the sopping criteria, c) led to incorrect choice of the damping parameter.

Setting scale alglib.minlmsetscale() solved this problem. I also recommend you to use less stringent stopping condition - epsg=1.0E-6 should be enough.

morbidlizard · **Joined:** Sat Jun 04, 2011 3:58 pm **Posts:** 5

I'm really sorry for being stupid and bothering you, but:
I tried to use minlmsetscale in such a way:

Code:

minlmsetscale(state, new double[]{10,10,10,10000});

right before invoking minlmoptimize(...)
but it didn't help - the solution is still not like the expected one. perhaps, i did sth wrong?
By the way, I'd like to know how to choose an appropriate vector for scaling(at least intuitively)
Thank you for your help

Sergey.Bochkanov · **Joined:** Fri May 07, 2010 7:06 am **Posts:** 927

1. Scale of the variable is its most typical value. You should know from problem formulation that one variable should be about 0.1-10 in magnitude while other can be as large as 100000.

2. Can you post exact (expected) answer? Your code at my computer returns

Quote:

X : 1,48064959286412, 0,805990620146673, 8,63291052919714, 254292,030484443
f(X) : 5,82076609134674E-11, 0,00818196564353979, 1,78589132451279E-05, -4,8129
8789711815E-08

and now I see that x[0] and x[1] are quite non-stable (they wander from 0 to their maximum values), and f[1] is significantly non-zero.

P.S. Don't worry to ask questions, I really appreciate feedback like yours.

morbidlizard · **Joined:** Sat Jun 04, 2011 3:58 pm **Posts:** 5

I got almost the same result - and f[1]'s not being zero is definitely wrong
thanks to mathcad, the expected answer is (approximately):

Code:

f[0] = 89.84
f[1] = 48.9
f[2] = 4.22
f[3] = 102374

morbidlizard · **Joined:** Sat Jun 04, 2011 3:58 pm **Posts:** 5

I've found something, that I consider to be a root of the problem.
I tried to calculate f for the approximate x (given in the message above) and have got an interesting result:

f(x) == 42.52 ; 0.038 ; almost 0 ; almost 0;

hey! f(x)[0] is not 0! but the solution is almost exact! it's not surprising that minlm* jumps over the correct solution and looks for it in wrong direction.
That led me to the idea of scaling f(x) too. The problem is (the new one), I don't know how to do it.
First of all I tried to divide f[i] by left part of corresponding equation, but it didn't help.
The source code now looks like

Code:

using System;
using System.Text;

namespace Project1
{
    class Program
    {
        //init
        static double d = 3;
        static double a = 1000;
        static double F_0 = 0.035;
        static double F_H = 549;
        static double F_d = 9;
        static double c_yd = 1.3;
        static double c_xd = 0.5;
        static double F_sh = 5;
        static double c_ysh = 0.6;
        static double c_xsh = 0.9;
        static double M_g = 1500;
        static double M_z = 400;
        static double hamma_M = 74000;
        static double V = 2.2;
        //
        static double ro = 1030;
        static double v = 1.3e-6;
        static double g = 9.80665;
        static double d_a = d / a;
        static double F_f = F_H * a / d;
        static double G_g = M_g * g * (hamma_M - ro * g) / hamma_M;
        static double G_z = M_z * g * (hamma_M - ro * g) / hamma_M;
        static double R_yd = c_yd * ro * V * V * F_d / 2;
        static double R_xd = c_xd * ro * V * V * F_d / 2;
        static double Q = c_ysh * ro * V * V * F_sh / 2;
        static double R_xsh = c_xsh * ro * V * V * F_sh / 2;
        
        //
        static void f(double[] x, double[] fi, object obj)
        {
            double L = x[0];// *x[0];
            double H = x[1];// * x[1];
            double A = x[2];// * x[2];
            double R_x = x[3];// * x[3];
            fi[0] = (R_x - (0.04 * (A + 110.0 * Math.Exp(-6.74 * F_0) * (Math.PI * H * L / (4.0 * F_f) - 0.02)) - 0.09) * (ro * V * V * F_H / 2.0));

            fi[1] = (Math.PI * H * L / (4.0 * F_f) - 1.0 + Math.Exp(0.0008 * (2.0 * R_yd + Q + 2.0 * G_g + G_z) / ((R_x + R_xsh) * F_0)));

            fi[2] = A + 110.0 * Math.Exp(-6.74 * F_0) * (Math.PI * H * L / (4.0 * F_f) - 0.02) - A * (1.0 - Math.Exp(-4.5 * (2.0 * R_yd + Q + 2.0 * G_g + G_z) / (R_x + R_xsh)));

            fi[3] = H / L - Math.Exp(1.0 - 2.0 * R_yd / (Q + 2.0 * G_g + G_z));
            return;
            //scaling
            fi[0] /= R_x;
            fi[1] /= H * L * Math.PI / (4 * F_f);
            fi[2] /= A * (1.0 - Math.Exp(-4.5 * (2.0 * R_yd + Q + 2.0 * G_g + G_z) / (R_x + R_xsh)));
            fi[3] /= H / L;
        }
        
        static double epsg = 1e-5;
        static double epsf = 0;
        static double epsx = 0;
        static int maxits = 0;
        //
        static void minlm(ref double[] x)
        {

            
            alglib.minlmstate state;
            alglib.minlmreport rep;
            alglib.minlmcreatev(4, x, 0.0001, out state);
            alglib.minlmsetcond(state, epsg, epsf, epsx, maxits);
            alglib.minlmsetscale(state, new double[] { 10, 10, 10, 10000 });
            alglib.minlmoptimize(state, f, null, null);
            alglib.minlmresults(state, out x, out rep);
            System.Console.WriteLine("{0}", rep.terminationtype); // EXPECTED: 4
        }
        static void Main()
        {




            double[] x = { 50, 50, 6, 10000 };

            double[] fi = { 0, 0, 0, 0 };

            minlm(ref x);
            f(x, fi, null);
            
            System.Console.WriteLine(" X : {0}, {1}, {2}, {3}", x[0], x[1], x[2], x[3]);
            System.Console.WriteLine(" f(X) : {0}, {1}, {2}, {3}", fi[0], fi[1], fi[2], fi[3]);
            System.Console.ReadLine();
        }
    }
}

I would appreciate any suggestions

Sergey.Bochkanov · **Joined:** Fri May 07, 2010 7:06 am **Posts:** 927

Quote:

hey! f(x)[0] is not 0! but the solution is almost exact!

Yes, I just noticed it. But it leads us to the question - are you sure that you've correctly converted your equation from Mathcad to C#? f0(x)=42 is too large to be considered "approximately zero" by Mathcad.

It can be the case that function your code actually calculates has minimum at different location than original one. And, unlike original problem, your incorrectly typed problem has no solution.

Can you focus on that possibility? and I will try to debug solver.

forum.alglib.net

Forum rules

nonlinear equations system

Who is online