Solving a nonlinear system of equations, with constraint of nonnegative solutions

I would like to solve a nonlinear system of 9 equations, with the constraint that the solutions must be nonnegative.

This problem arises in the context of numerically solving for an Endemic Equilibrium, an equilibrium where both diseases in a population are present.

I have tried with the code below - but I have a feeling something is wrong. Even though I get a result, it is not biologically realistic. I should have positive values for all nine solutions, but am getting the following (see result below).

Can someone help me?

Input:

function test()

xo = [191564208 131533276 2405659 1805024 1000000 1000000 500000 500000]

Lambda = 531062; 
mu = (1/70)/365;
mu_A = 0.25/365;
mu_T = 0.2/365;
beta = 0.187/365; 
tau = 4/365;  
eta_1 = 0;
eta_2 = 1;
lambda_T = 0.1;
rho_1 = 1/60;
rho_2 = (rho_1)./(270.*rho_1-1);
gamma = 1e-3;

      options = optimset('MaxFunEvals',1E4, ...
          'MaxIter', 1E4, ...
          'TolFun', 1E-32, ...
          'TolX', 1E-32, ...
          'TolCon', 1E-32);

      x = fmincon(@(X) Ftest(X), xo, [], [], [], [], ...
          [0 0 0 0 0 0 0 0], [], @(X) xcon(X),  options)

      function F = Ftest(x)
          F = norm(Research(x))
      end

function [c,ceq] = xcon(x)
        c = []
        ceq = [Lambda-mu.*x(1)-beta.*(x(3)+x(4)+x(5)+x(6)).*(x(1)./(x(1)+x(2)+x(3)+x(4)+x(5)+x(6)+x(7)+x(8)))-tau.*(x(2)+x(4)).*(x(1)./(x(1)+x(2)+x(3)+x(4)+x(5)+x(6)+x(7)+x(8)));
            tau.*(x(2)+x(4)).*(x(1)./(x(1)+x(2)+x(3)+x(4)+x(5)+x(6)+x(7)+x(8)))-beta.*(x(3)+x(4)+x(5)+x(6)).*(x(2)./(x(1)+x(2)+x(3)+x(4)+x(5)+x(6)+x(7)+x(8)))-(mu+mu_T).*x(2);
            beta.*(x(3)+x(4)+x(5)+x(6)).*(x(1)./(x(1)+x(2)+x(3)+x(4)+x(5)+x(6)+x(7)+x(8)))-tau.*(x(2)+x(4)).*(x(3)./(x(1)+x(2)+x(3)+x(4)+x(5)+x(6)+x(7)+x(8)))-(mu+mu_A).*x(3);
            beta.*(x(3)+x(4)+x(5)+x(6)).*(x(2)./(x(1)+x(2)+x(3)+x(4)+x(5)+x(6)+x(7)+x(8)))+tau.*(x(2)+x(4)).*(x(3)./(x(1)+x(2)+x(3)+x(4)+x(5)+x(6)+x(7)+x(8)))-(mu+mu_A+mu_T+lambda_T).*x(4);
            lambda_T.*x(4)-(mu+mu_A+rho_1+eta_1).*x(5);
            rho_1.*x(5)-(mu+mu_A+rho_2+eta_2).*x(6);
            eta_1.*x(5)-(mu+rho_1+gamma).*x(7);
            eta_2.*x(6)-(mu+rho_2+gamma.*(rho_1)./(rho_1+rho_2)).*x(8)+(rho_1).*x(7)];

end

      function F = Research(x)
          F = [Lambda-mu.*x(1)-beta.*(x(3)+x(4)+x(5)+x(6)).*(x(1)./(x(1)+x(2)+x(3)+x(4)+x(5)+x(6)+x(7)+x(8)))-tau.*(x(2)+x(4)).*(x(1)./(x(1)+x(2)+x(3)+x(4)+x(5)+x(6)+x(7)+x(8)));
              tau.*(x(2)+x(4)).*(x(1)./(x(1)+x(2)+x(3)+x(4)+x(5)+x(6)+x(7)+x(8)))-beta.*(x(3)+x(4)+x(5)+x(6)).*(x(2)./(x(1)+x(2)+x(3)+x(4)+x(5)+x(6)+x(7)+x(8)))-(mu+mu_T).*x(2);
              beta.*(x(3)+x(4)+x(5)+x(6)).*(x(1)./(x(1)+x(2)+x(3)+x(4)+x(5)+x(6)+x(7)+x(8)))-tau.*(x(2)+x(4)).*(x(3)./(x(1)+x(2)+x(3)+x(4)+x(5)+x(6)+x(7)+x(8)))-(mu+mu_A).*x(3);
              beta.*(x(3)+x(4)+x(5)+x(6)).*(x(2)./(x(1)+x(2)+x(3)+x(4)+x(5)+x(6)+x(7)+x(8)))+tau.*(x(2)+x(4)).*(x(3)./(x(1)+x(2)+x(3)+x(4)+x(5)+x(6)+x(7)+x(8)))-(mu+mu_A+mu_T+lambda_T).*x(4);
              lambda_T.*x(4)-(mu+mu_A+rho_1+eta_1).*x(5);
              rho_1.*x(5)-(mu+mu_A+rho_2+eta_2).*x(6);
              eta_1.*x(5)-(mu+rho_1+gamma).*x(7);
              eta_2.*x(6)-(mu+rho_2+gamma.*(rho_1)./(rho_1+rho_2)).*x(8)+(rho_1).*x(7)];

      end

end

Output:

Local minimum possible. Constraints satisfied.

fmincon stopped because the size of the current search direction is less than
twice the selected value of the step size tolerance and constraints are 
satisfied to within the selected value of the constraint tolerance.

<stopping criteria details>

Active inequalities (to within options.TolCon = 1e-32):
  lower      upper     ineqlin   ineqnonlin
    3                                 
    4                                 
    5                                 
    6                                 
    7                                 
    8                                 

x =

     1.0e+08 *

    Columns 1 through 4

     0.510099026315789   9.011749464912279  -0.000000000000000  -0.000000000000000

    Columns 5 through 8

    -0.000000000000000  -0.000000000000000                   0  -0.000000000000000

[Torsten]

1. I only see eight equations.

2. I don't see the objective function Ftest.

3. You can easily test whether the result is at least numerically correct by evaluating your equality constraints for the obtained solution..

Best wishes

Torsten.

[Torsten]

Sorry, now I see Ftest.

But since you only want the constraints to be satisfied, you can also use FTest = @(x)0.

Best wishes

Torsten.

Answer 1

[Matt J]

As the exit message tells you, your constraints are satisfied within the constraint tolerance that you supplied, TolCon=1e-32.

If you don't truly require the bounds to be satisfied exactly, then fmincon was successful and you have your solution. However, using fmincon in this way is massive overkill. You could just as well have used lsqnonlin.

If you do require the bounds to be exactly satisfied, you should use fmincon's interior-point or sqp algorithm with their default settings. Note however, that lb,ub bounds are the only kind of constraints that can be satisfied exactly. Any other constraint can only be guaranteed to within TolCon and, in the case of nonlinear constraints, only when you provide a sufficiently good initial guess, xo.

Also, similar to Torsten comments, it should be unnecessary to have both nonlinear constraints and a least squares cost function that are identical to each other. However, I think it makes more sense to get rid of xcon. If a solution to your equations exists, it should be enough to minimize norm(Research(x)). Since nonlinear constraints can't be satisfied exactly, there is no additional benefit to including xcon. It just adds computational burden.