Nonlinear Programming
The goal of nonlinear programming is to optimize a possibly nonlinear function subject to linear or nonlinear constraints. In its most general form, a nonlinear program is an optimization problem
Minimize | f(x) |
|---|---|
subject to | lhsi ≤ Aix ≤ rhsi |
lhsj ≤ gj(x) ≤ rhsj | |
li ≤ xi ≤ ui |
The function f is a nonlinear function of the xi. The matrix A contains the coefficients of the linear constraints. The functions gj represent nonlinear constraint functions.
Most nonlinear programs are relatively small, but some have hundreds or thousands of variables and constraints.
Variables and Constraints
The DecisionVariable class is used to represent variables in a nonlinear program. Each variable has a unique name, available through the Name property. The variables are the unknowns in the optimization problem. The objective function is a nonlinear function of the variables.
The LowerBound and UpperBound properties define constraints on the values of a variable. Standard variables have no lower or upper bound. When a variable is not bounded from above, its upper bound is Double.PositiveInfinity. When a variable is not bounded from below, its lower bound is Double.NegativeInfinity.
The Constraint class is the abstract base class used to represent constraints in a nonlinear program. Every constraint must have a unique name, which can be accessed through the Name property. The LowerBound and UpperBound properties define these boundaries. When a constraint is not bounded from above, its upper bound is Double.PositiveInfinity. When a constraint is not bounded from below, its lower bound is Double.NegativeInfinity.
A linear constraint is a linear combination of variables that is required to lie between specified boundaries. Linear constraints are represented by the LinearConstraint class.
A nonlinear constraint is a (possibly) nonlinear function that is required to lie between specified boundaries. Nonlinear constraints are represented by the NonlinearConstraint class.
Variables and constraints can be accessed through the nonlinear program's Variables and Constraints collections. These collections are indexed by position and by name.
Defining nonlinear programs
A nonlinear program is represented by the NonlinearProgram class. There are two ways to construct NonlinearProgram objects.
Defining nonlinear programs using constructors
The quickest way to define a nonlinear program is by passing the vectors and matrices that appear in the definition above to a constructor. There are three constructors that can be used in this way.
The first constructor defines a nonlinear program with only linear constraints in standard form. This is a nonlinear program where variables are constrained to be positive, and all constraints are inequalities with an upper bound. It takes four arguments. The first argument is a Func<T, TResult> that maps a vector to a real number and represents the objective function. The second argument is a Func<T1, T2, TResult> that represents the gradient of the objective function. The first argument is a vector containing the x values. The second argument is a vector that is to hold the result.
The third argument is a Matrix<T> that contains the coefficients of the inequalities. It corresponds to the matrix A in the definition. The fourth argument is a Vector<T> containing the right-hand-sides of the inequalities. The code below creates the following a nonlinear program:
Minimize | x2 + 4y2 - 32y + 64 |
|---|---|
s.t. | x + y ≤ 7 |
-x + 2y ≤ 4 | |
x ≥ 0 | |
y ≥ 0 |
Func<Vector<double>, double> f = x =>
x[0] * x[0] + 4 * x[1] * x[1] - 32 * x[1] + 64;
Func<Vector<double>, Vector<double>, Vector<double>>
g = (x, y) =>
{
y[0] = 2 * x[0];
y[1] = 8 * x[1] - 32;
return y;
};
var A = Matrix.CopyFrom(new[,] { { 1.0, 1.0 }, { -1.0, 2.0 } });
var b = Vector.Create(7.0, 4.0);
var nlp1 = new NonlinearProgram(f, g, A, b, 0);The third constructor is the most flexible. It has seven parameters in all. The first two are once again the objective function and its gradient. The third argument is the coefficient matrix. The fourth and fifth parameters are vectors that contain the lower and upper bounds of the constraints. Each type of constraint can be expressed using a lower and upper bound, as shown in the table below.
Constraint | Lower bound | Upper bound |
|---|---|---|
c ≤ rhs | -infinity | rhs |
c ≥ rhs | rhs | +infinity |
c = rhs | rhs | rhs |
lhs ≤ c ≤ rhs | lhs | rhs |
The sixth and seventh parameters are vectors that contain the lower and upper bound for the variables. The different types of constraints on variables can be expressed the same way as for constraints. The following code sample builds the same nonlinear program:
Func<Vector<double>, double> f = x =>
x[0] * x[0] + 4 * x[1] * x[1] - 32 * x[1] + 64;
Func<Vector<double>, Vector<double>, Vector<double>>
g = (x, y) =>
{
y[0] = 2 * x[0];
y[1] = 8 * x[1] - 32;
return y;
};
var A = Matrix.CopyFrom(new[,] { { 1.0, 1.0 }, { -1.0, 2.0 } });
var lc = Vector.CreateConstant(2, double.NegativeInfinity);
var uc = Vector.Create(7.0, 4.0);
var lv = Vector.CreateConstant(2, 0.0);
var uv = Vector.CreateConstant(2, double.PositiveInfinity);
var nlp2 = new NonlinearProgram(f, g, A, lc, uc, lv, uv);Building nonlinear programs from scratch
The NonlinearProgram class also has two constructors that creates an empty nonlinear program. If you don't pass any arguments, the model is completely empty and you have to explicitly add variables to it. You can then call the nonlinear program's AddVariable methods to define variables. You may also pass in the number of variables, and a set of standard variables will be created automatically.
The AddVariable method has two overloads. The first has one parameter: the name of the variable. The second overload has two additional parameters that specify the lower and upper bound of the variable. If no values are specified, a lower bound of negative infinity and an upper bound of infinity are assumed.
The AddLinearConstraint method has four overloads. The first takes four parameters. The first is the name of the constraint. The second is a Vector<T> that contains the coefficients of the variables. The third parameter is a ConstraintType value that specifies the type of constraint. The possible values are listed in the table below. The last parameter is the right-hand side of the constraint.
Value | Description |
|---|---|
Equal | c = rhs |
GreaterThanOrEqual | c ≥ rhs |
LessThanOrEqual | c ≤ rhs |
The second overload also takes four parameters. The first two are the same as before. The third and fourth parameters are the lower and upper bound of the constraint. The different types of constraints can be specified according to the examples in Table 1.
The third and fourth overloads are similar to the first two. The only difference is that the coefficients are specified using a Double array instead of a vector.
The following example creates a nonlinear program that is identical to the last example:
var nlp3 = new NonlinearProgram();
nlp3.AddVariable("x1", 0.0, double.PositiveInfinity);
nlp3.AddVariable("x2", 0.0, double.PositiveInfinity);
nlp3.ObjectiveFunction =
x => x[0] * x[0] + 4 * x[1] * x[1] - 32 * x[1] + 64;
nlp3.ObjectiveGradient =
(x, y) =>
{
y[0] = 2 * x[0];
y[1] = 8 * x[1] - 32;
return y;
};
nlp3.AddLinearConstraint("c1", new[] { 1.0, 1.0 },
ConstraintType.LessThanOrEqual, 7.0);
nlp3.AddLinearConstraint("c1", new[] { -1.0, 2.0 },
ConstraintType.LessThanOrEqual, 4.0);Although not strictly necessary, it is advisable to create the variables first, followed by the constraints. If the vector containing the coefficients for a constraint is longer than the number of variables that have been defined, additional variables are added automatically. The kth variable is given the name x<k>.
Whenever possible, constraints on variables should be preferred over explicit nonlinear program constraints. The time to solve a nonlinear program depends most heavily on the number of constraints.
Using Automatic Differentiation
The gradient of the objective function and nonlinear constraints are needed during the solution process. When using .NET 4.0, it is possible to have these computed automatically using a symbolic representation of the objective function or constraint.
To do so, create the nonlinear program, and then set its SymbolicObjectiveFunction to a lambda expression that represents the objective function. The lambda expression should take a Vector<T> as its input. To define nonlinear constraints, call one of the overloads of AddSymbolicConstraint.
Choose a Method and Supply Second Derivatives
Create a NonlinearProgramOptions to configure the method, tolerances and step policies, then pass it to Solve(NonlinearProgramOptions). An overload also accepts parallel and cancellation settings, analogously to linear-program solves. Each solve copies the options before any derivative or objective callback. Changing the supplied object affects a subsequent solve, not the active run. A parameterless solve initially uses defaults; after passing options, it takes a fresh copy of the most recently supplied object. Derivative callbacks, structure and the initial guess remain on the model. Iteration and evaluation limits remain on its Solver.
The default Method is SequentialQuadraticProgramming. Keep it as the starting choice for a general dense nonlinear program. ReducedSequentialQuadraticProgramming is an alternative that uses curvature along directions allowed by the active constraints. This can be useful when the Hessian, the matrix of second derivatives, is indefinite even though curvature along feasible directions is positive.
Set Method to Automatic to select once at initialization. A supplied sparse constraint Jacobian selects sparse SQP; other generally constrained models use dense SQP. Bounds-only models use Newton-CG when a Hessian is supplied, otherwise bounded or unbounded limited-memory BFGS. Explicit methods take precedence. Automatic selection updates the supplied options' Method to the concrete choice before objective and Hessian callbacks. Reusing those options retains the choice; set Method to Automatic again to request new selection. A method change made by the caller during the initial constraint-Jacobian callback is preserved for the next solve, while the active solve continues using its own snapshot. The default remains SequentialQuadraticProgramming.
Both methods use Hessian when you supply it. This function includes second derivatives of the objective and the constraints, weighted by the multipliers from the current step. The API reference gives the required signs and ordering. When you omit it, the solver builds a BFGS approximation from changes in first derivatives. Gradients are approximated numerically when you do not supply them.
CurvatureDampingFactor controls how much positive curvature is retained when a gradient change indicates weak or negative curvature. This uses Powell damping: it blends the gradient change with the current Hessian times the step before the BFGS update. It does not add a diagonal shift. Set a finite value strictly between zero and one before solving. A larger value causes more damping. The default is null, which preserves 0.1 for both dense SQP methods and 0.2 for sparse compact BFGS. Supplied Hessians and inverse-history preconditioning keep their own policies.
The following example supplies first and second derivatives and explicitly selects the reduced method. Its constraint is linear, so only the objective contributes second derivatives.
NonlinearProgram model = new NonlinearProgram(2);
model.InitialGuess = Vector.Create(.5, .5);
model.ObjectiveFunction = x => (x[0]-1)*(x[0]-1) + (x[1]-2)*(x[1]-2);
model.ObjectiveGradient = (x,g) =>
{
g[0] = 2*(x[0]-1); g[1] = 2*(x[1]-2);
return g;
};
model.AddLinearConstraint("total", Vector.Create(1.0, 1), ConstraintType.LessThanOrEqual, 2);
// Optional: supply second derivatives when they are available.
// Without this function, the solver builds a BFGS approximation.
model.Hessian = (x,weight,multipliers,previous) =>
{
Matrix<double> h = previous as Matrix<double> ?? Matrix.CreateSymmetric<double>(2);
h.Clear();
h[0,0] = 2*weight; h[1,1] = 2*weight;
return h;
};
NonlinearProgramOptions options = new NonlinearProgramOptions
{ Method = NonlinearProgramMethod.ReducedSequentialQuadraticProgramming, OptimalityTolerance = 1e-6 };
Vector<double> solution = model.Solve(options);
if (model.SolutionReport.Status != AlgorithmStatus.Converged)
throw new InvalidOperationException("The nonlinear solve did not converge.");
// The solution is (0.5, 1.5).Set OptimalityTolerance to adjust the final checks on constraints and first-order optimality. The default is 1e-6. Constraint violations are checked in the units of your model; derivative checks use relative scales. A smaller tolerance can require more iterations or cause failure if the steps cannot achieve enough accuracy. A small step alone does not establish convergence, and convergence does not prove global optimality for a nonconvex problem.
UseQuadraticProgramWarmStart lets the reduced method reuse active constraints between successive QP steps in one nonlinear solve. It is off by default. It can save work when successive steps have similar active constraints, but preparing and checking the starting point also takes time. Measure it on representative problems before enabling it. Each new nonlinear solve resets this state. Changes to the supplied method, tolerance or warm-start options affect the next solve; the current solve retains its initial copy.
Methods for Large Problems
Dense SQP stores matrices whose size grows with the square of the number of variables. For large problems, choose a method that uses limited history, matrix products or sparse entries. The best choice depends on the constraints and the derivatives you can provide.
LimitedMemoryBfgs is for problems without bounds or general constraints. BoundedLimitedMemoryBfgs also handles variable bounds. These methods use changes in gradients to build a small history of useful directions. They do not require second derivatives. Set the history size through HistorySize; the default is 10. They are a useful starting choice when only first derivatives are available.
NewtonConjugateGradient supports variable bounds but no general constraints. It needs products of the Hessian with a direction, rather than the full matrix. Supply Hessian when available. The callback returns a linear operator and receives an objective weight and multipliers. For Newton's bounds-only problems the multiplier vector is empty. The objective weight is 1 for minimization and -1 for maximization, consistently with SQP. The operator supplies Hessian-vector multiplication. The following example has 1,000 variables and supplies products without creating a square matrix.
const int n = 1000;
NonlinearProgram model = new NonlinearProgram(n);
model.InitialGuess = Vector.Zeros<double>(n);
for (int i = 0; i < n; i++) { model.Variables[i].LowerBound = 0; model.Variables[i].UpperBound = 2; }
// Minimize sum((x[i]-1)^2)/2. The solution has every component equal to 1.
model.ObjectiveFunction = x =>
{ double value = 0; for (int i = 0; i < n; i++) value += .5*(x[i]-1)*(x[i]-1); return value; };
model.ObjectiveGradient = (x,g) => { for (int i = 0; i < n; i++) g[i] = x[i]-1; return g; };
// H is the identity: multiplying by H just copies the direction.
model.Hessian = (x, weight, multipliers, previous) => previous ??
LinearOperator<double>.Create(n, n, (operation, direction, result) =>
{
result = result ?? Vector.Zeros<double>(n);
for (int i = 0; i < n; i++) result[i] = weight * direction[i];
return result;
});
model.HessianDiagonal = (x,r) => { for (int i = 0; i < n; i++) r[i] = 1; };
NonlinearProgramOptions options = new NonlinearProgramOptions
{ Method = NonlinearProgramMethod.NewtonConjugateGradient };
Vector<double> solution = model.Solve(options);
if (model.SolutionReport.Status != AlgorithmStatus.Converged) throw new InvalidOperationException("Newton did not converge.");Newton can also use sparse Hessian entries, or approximate products by differencing gradients. If no gradients are supplied, it approximates those from objective values too. This avoids dense storage but can require many function calls. Numerical products can be less reliable when function values or gradients contain noise; compare with the limited-memory method. Bound-only numerical differences stay inside the variable bounds.
A preconditioner makes Newton's inner linear solve easier. The automatic choice uses incomplete Cholesky for sparse entries, a supplied diagonal when available, or limited-memory history. Incomplete Cholesky builds an inexpensive approximate factor of a sparse Hessian. ModifiedIncompleteCholesky is an explicit alternative for some grid problems. Unbounded Newton problems can optionally use a backtracking line search through UseLineSearch.
An optional HessianDiagonal callback supplies the objective Hessian diagonal for Newton preconditioning. It receives the point and a destination vector; overwrite every component, including entries that become zero. Use the original objective sign, including when maximizing, without the weight or constraint terms received by the Hessian callback. An explicit Hessian matrix can supply its diagonal when this callback is omitted.
For an objective written as a sum of squared residuals, you can explicitly supply Gauss-Newton curvature, 2 JTJ, through a Hessian matrix or multiplication operator. Here J contains derivatives of the residuals. This approximation omits their second derivatives; the objective and gradient callbacks still evaluate the original function. For some very flat problems, this curvature combined with UseLineSearch converges more reliably than exact-Hessian trust-region steps. Supply this structure explicitly and compare it on your models. A small objective or gradient does not guarantee accurate coordinates on a nearly flat objective.
SparseSequentialQuadraticProgramming supports general constraints. Supply ConstraintJacobianFunction, returning a sparse matrix, with an optional Hessian function returning an explicit matrix. The Jacobian contains first derivatives of the constraints. Its row indexes follow the Constraints collection; column indexes follow the Variables collection. The first call receives a null reusable matrix. Later calls can update that matrix or return a replacement. Overwrite entries that change, including entries becoming zero. Stored structure must include entries that can become nonzero later, even when they are zero at the initial point. Numerical nonzero enumeration cannot establish this pattern. For an initially zero structural entry, create the stored position first and protect the matrix with MutableValues mutability before writing zero. This preserves the position while allowing subsequent value updates. Ordinary sparse constructors and mutable-structure updates may remove zeros.
Hessian values describe the Lagrangian using the same weights and signs as the dense callback. Return a symmetric matrix. The initial stored structure is captured, including stored zeros; later structural growth is rejected. Without a Hessian function, an optional HessianPattern matrix template permits numerical recovery from first derivatives. Without a Hessian pattern, the method uses compact BFGS. Its damping uses the same CurvatureDampingFactor setting described above. It does not create a dense Jacobian as a fallback. These callbacks must work at trial points, including points outside variable bounds for this infeasible-start method.
Hessian functions receive a null reusable operator on the first call, and may update or replace it on subsequent calls. The operator must remain valid for all products until the next callback; copy any input data retained by an operator. Newton supports multiplication-only operators. Direct SQP methods require explicit matrix entries and reject opaque operators rather than reconstructing dense matrices. Sparse IC/MIC preconditioning also requires explicit sparse entries.
NonlinearProgram model = new NonlinearProgram(3);
model.InitialGuess = Vector.Zeros<double>(3);
model.ObjectiveFunction = x => .5*((x[0]-1)*(x[0]-1)+(x[1]-2)*(x[1]-2)+(x[2]-3)*(x[2]-3));
model.ObjectiveGradient = (x,g) => { for (int i = 0; i < 3; i++) g[i] = x[i]-(i+1); return g; };
model.AddNonlinearConstraint("first",x => x[0]+x[1],ConstraintType.Equal,2);
model.AddNonlinearConstraint("second",x => x[1]+x[2],ConstraintType.Equal,4);
model.ConstraintJacobianFunction = (x, previous) =>
{
Matrix<double> jacobian = previous ?? Matrix.CreateSparse<double>(2, 3);
jacobian[0, 0] = 1; jacobian[0, 1] = 1;
jacobian[1, 1] = 1; jacobian[1, 2] = 1;
return jacobian;
};
model.Hessian = (x, weight, multipliers, previous) =>
{
Matrix<double> hessian = previous as Matrix<double> ?? Matrix.CreateSparseSymmetric<double>(3);
// These constraints are linear, so their second derivatives are zero.
for (int i = 0; i < 3; i++) hessian[i, i] = weight;
return hessian;
};
NonlinearProgramOptions options = new NonlinearProgramOptions
{ Method = NonlinearProgramMethod.Automatic };
// The initial sparse Jacobian selects sparse SQP and updates options.Method.
Vector<double> solution = model.Solve(options);
if (model.SolutionReport.Status != AlgorithmStatus.Converged) throw new InvalidOperationException("Sparse optimization did not converge.");
// The solution is (2/3, 4/3, 8/3).Sparse optimization uses adaptive factorization by default. It can try a short feasibility recovery phase and constraint-curvature step corrections when ordinary steps fail. UseFeasibilityRestoration and UseSecondOrderCorrection let you disable either recovery technique. Sparse storage still depends on fill created during factorization; it is not guaranteed to grow linearly for every constraint graph.
Set Solver.MaxIterations and Solver.MaxEvaluations to limit work. Large-scale methods count objective calls made during numerical derivative probes as well as ordinary trials. Hessian-vector products alone do not count as objective evaluations. Cancellation through Solver.ParallelOptions raises OperationCanceledException and preserves the available point for inspection. All nonlinear methods honor MinIterations and ConvergenceTests at iteration checkpoints. Additional tests can delay successful termination, but cannot replace the original optimality checks. Their enabled flags, collection quantifier, pre-iteration behavior and required successive passes apply. Successive-pass counts reset for each solve. An already stationary point is retained while satisfying additional iteration or test requirements; the maximum iteration and evaluation limits still apply. Custom tests can inspect the current solver report, and their failure statuses or exceptions are preserved. Every converged result must pass fresh final checks at OptimalityTolerance. A failed sparse solve does not establish infeasibility. Convergence on a nonconvex objective does not prove either a global optimum or a local minimum. Solve-local histories and derivative caches reset on each call; changing methods preserves the Solver object and its limits.
Solving nonlinear programs
Once a nonlinear program has been defined, solving it is very straightforward. The Solve method does all the work. Sequential quadratic programming (SQP) solves a sequence of quadratic approximations to the problem. It checks each proposed step against the original objective and constraints before accepting it. The two dense SQP choices are intended for small and moderate problems. The explicit large-problem methods described above use less storage.
The Solve method returns a Vector<T> containing the optimal solution. You should inspect the Status property to make sure an optimal solution was indeed found. This property is of type OptimizationModelStatus and can take on the following values:
Value | Description |
|---|---|
Unknown | Nothing is known about the solution. This is the status returned before the SolveModel method is called. |
Optimal | The solver converged and the final solution checks passed. For a nonconvex problem, this does not prove global optimality. |
Infeasible | The nonlinear program does not have a solution because some of the constraints conflict with each other. |
Unbounded | The nonlinear program does not have a finite solution. The cost function can be made arbitrarily small. |
Finally, the OptimalValue property returns the value of the cost function at the solution.
The code below solves the nonlinear program defined earlier:
var solution = nlp3.Solve();
Console.WriteLine(solution);