Nonlinear Systems with Constraints

Solve Equations with Inequality Constraints

fsolve solves a system of nonlinear equations. However, it does not allow you to include any constraints, even bound constraints. So how can you solve a system of nonlinear equations when you have constraints?

A solution that satisfies your constraints is not guaranteed to exist. In fact, the problem might not have any solution, even one that does not satisfy your constraints. However, techniques exist to help you search for solutions that satisfy your constraints.

To illustrate the techniques, consider how to solve the equations

where the components of $$x$$ must be nonnegative. The equations have four solutions:

Only one solution satisfies the constraints, namely $$x = (10,20)$$.

The fbnd helper function at the end of this example calculates $$F(x)$$ numerically.

Use Different Start Points

Generally, a system of $$N$$ equations in $$N$$ variables has isolated solutions, meaning each solution has no nearby neighbors that are also solutions. So, one way to search for a solution that satisfies some constraints is to generate a number of initial points x0, and then run fsolve starting at each x0.

For this example, to look for a solution to the equation system $$F(x)=0$$, take 10 random points that are normally distributed with mean 0 and standard deviation 100.

rng default % For reproducibility
N = 10; % Try 10 random start points
pts = 100*randn(N,2); % Initial points are rows in pts
soln = zeros(N,2); % Allocate solution
opts = optimoptions('fsolve','Display','off');
for k = 1:N
    soln(k,:) = fsolve(@fbnd,pts(k,:),opts); % Find solutions
end

List solutions that satisfy the constraints.

idx = soln(:,1) >= 0 & soln(:,2) >= 0;
disp(soln(idx,:))

   10.0000   20.0000
   10.0000   20.0000
   10.0000   20.0000
   10.0000   20.0000
   10.0000   20.0000

Use Different Algorithms

fsolve has three algorithms. Each can lead to different solutions.

For this example, take x0 = [1,9] and examine the solution each algorithm returns.

x0 = [1,9];
opts = optimoptions(@fsolve,'Display','off',...
    'Algorithm','trust-region-dogleg');
x1 = fsolve(@fbnd,x0,opts)

x1 = 1×2

   -1.0000   -2.0000

opts.Algorithm = 'trust-region';
x2 = fsolve(@fbnd,x0,opts)

x2 = 1×2

   -1.0000   20.0000

opts.Algorithm = 'levenberg-marquardt';
x3 = fsolve(@fbnd,x0,opts)

x3 = 1×2

    0.9523    8.9941

Here, all three algorithms find different solutions for the same initial point. None satisfy the constraints. The reported "solution" x3 is not even a solution, but is simply a locally stationary point.

Use lsqnonlin with Bounds

lsqnonlin tries to minimize the sum of squares of the components in a vector function $$F(x)$$. Therefore, it attempts to solve the equation $$F(x) = 0$$. Also, lsqnonlin accepts bound constraints.

Formulate the example problem for lsqnonlin and solve it.

lb = [0,0];
rng default
x0 = 100*randn(2,1);
[x,res] = lsqnonlin(@fbnd,x0,lb)

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

x = 2×1

   10.0000
   20.0000

res = 
2.4783e-25

In this case, lsqnonlin converges to the solution satisfying the constraints. You can use lsqnonlin with the Global Optimization Toolbox MultiStart solver to search over many initial points automatically. See MultiStart Using lsqcurvefit or lsqnonlin (Global Optimization Toolbox).

Set Equations and Inequalities as fmincon Constraints

You can reformulate the problem and use fmincon as follows:

The fminconstr helper function at the end of this example implements the nonlinear constraints. Solve the constrained problem.

lb = [0,0]; % Lower bound constraint
rng default % Reproducible initial point
x0 = 100*randn(2,1);
opts = optimoptions(@fmincon,'Algorithm','interior-point','Display','off');
x = fmincon(@(x)0,x0,[],[],[],[],lb,[],@fminconstr,opts)

x = 2×1

   10.0000
   20.0000

In this case, fmincon solves the problem from the start point.

Helper Functions

This code creates the fbnd helper function.

function F = fbnd(x)

F(1) = (x(1)+1)*(10-x(1))*(1+x(2)^2)/(1+x(2)^2+x(2));
F(2) = (x(2)+2)*(20-x(2))*(1+x(1)^2)/(1+x(1)^2+x(1));
end

This code creates the fminconstr helper function.

function [c,ceq] = fminconstr(x)

c = []; % No nonlinear inequality
ceq = fbnd(x); % fsolve objective is fmincon nonlinear equality constraints
end