mpcqpsolver

(To be removed) Solve a quadratic programming problem using the KWIK algorithm

mpcqpsolver will be removed in a future release. Use mpcActiveSetSolver instead. For more information, see Version History.

Syntax

[x,status] = mpcqpsolver(Linv,f,A,b,Aeq,beq,iA0,options)

[x,status,iA,lambda] = mpcqpsolver(Linv,f,A,b,Aeq,beq,iA0,options)

Description

[x,status] = mpcqpsolver(Linv,f,A,b,Aeq,beq,iA0,options) finds an optimal solution, x, to a quadratic programming problem by minimizing the objective function:

$$J=\frac{1}{2}{x}^{⊺}Hx+f^{⊺}x$$

subject to inequality constraints $$Ax\ge b$$, and equality constraints $$A_{eq}x=b_{eq}$$. status indicates the validity of x.

[x,status,iA,lambda] = mpcqpsolver(Linv,f,A,b,Aeq,beq,iA0,options) also returns the active inequalities, iA, at the solution, and the Lagrange multipliers, lambda, for the solution.

Examples

Solve Quadratic Programming Problem Using Active-Set Solver

Find the values of x that minimize

$$f\left(x\right)=0.5x_1^2+x_2^2-x_1{x}_2-2x_1-6x_2,$$

subject to the constraints

$$\begin{array}{l}{x}_1\ge 0\\ x_2\ge 0\\ x_1+x_2\le 2\\ -x_1+2x_2\le 2\\ 2x_1+x_2\le 3.\end{array}$$

Specify the Hessian and linear multiplier vector for the objective function.

H = [1 -1; -1 2];
f = [-2; -6];

Specify the inequality constraint parameters.

A = [1 0; 0 1; -1 -1; 1 -2; -2 -1];
b = [0; 0; -2; -2; -3];

Define Aeq and beq to indicate that there are no equality constraints.

Aeq = [];
beq = zeros(0,1);

Find the lower-triangular Cholesky decomposition of H.

[L,p] = chol(H,'lower');
Linv = inv(L);

It is good practice to verify that H is positive definite by checking if p = 0.

p

p = 0

Create a default option set for mpcActiveSetSolver.

opt = mpcqpsolverOptions;

To cold start the solver, define all inequality constraints as inactive.

iA0 = false(size(b));

Solve the QP problem.

[x,status] = mpcqpsolver(Linv,f,A,b,Aeq,beq,iA0,opt);

Examine the solution, x.

x

x = 2×1

    0.6667
    1.3333

Check Active Inequality Constraints for QP Solution

Find the values of x that minimize

$$f\left(x\right)=3x_1^2+0.5x_2^2-2x_1{x}_2-3x_1+4x_2,$$

subject to the constraints

$$\begin{array}{l}{x}_1\ge 0\\ x_1+x_2\le 5\\ x_1+2x_2\le 7.\end{array}$$

Specify the Hessian and linear multiplier vector for the objective function.

H = [6 -2; -2 1];
f = [-3; 4];

Specify the inequality constraint parameters.

A = [1 0; -1 -1; -1 -2];
b = [0; -5; -7];

Define Aeq and beq to indicate that there are no equality constraints.

Aeq = [];
beq = zeros(0,1);

Find the lower-triangular Cholesky decomposition of H.

[L,p] = chol(H,'lower');
Linv = inv(L);

Verify that H is positive definite by checking if p = 0.

p

p = 0

Create a default option set for mpcqpsolver.

opt = mpcqpsolverOptions;

To cold start the solver, define all inequality constraints as inactive.

iA0 = false(size(b));

Solve the QP problem.

[x,status,iA,lambda] = mpcqpsolver(Linv,f,A,b,Aeq,beq,iA0,opt);

Check the active inequality constraints. An active inequality constraint is at equality for the optimal solution.

iA

iA = 3x1 logical array

   1
   0
   0

There is a single active inequality constraint.

View the Lagrange multiplier for this constraint.

lambda.ineqlin(1)

ans = 5.0000

Input Arguments

Linv — Inverse of lower-triangular Cholesky decomposition of Hessian matrix

n-by-n matrix

Inverse of lower-triangular Cholesky decomposition of Hessian matrix, specified as an n-by-n matrix, where n > 0 is the number of optimization variables. For a given Hessian matrix, H, Linv can be computed as follows:

[L,p] = chol(H,'lower');
Linv = inv(L);

H is an n-by-n matrix, which must be symmetric and positive definite. If p = 0, then H is positive definite.

Note

The KWIK algorithm requires the computation of Linv instead of using H directly, as in the quadprog (Optimization Toolbox) command.

f — Multiplier of objective function linear term

column vector

Multiplier of objective function linear term, specified as a column vector of length n.

A — Linear inequality constraint coefficients

m-by-n matrix | []

Linear inequality constraint coefficients, specified as an m-by-n matrix, where m is the number of inequality constraints.

If your problem has no inequality constraints, use [].

b — Right-hand side of inequality constraints

column vector of length m

Right-hand side of inequality constraints, specified as a column vector of length m.

If your problem has no inequality constraints, use zeros(0,1).

Aeq — Linear equality constraint coefficients

q-by-n matrix | []

Linear equality constraint coefficients, specified as a q-by-n matrix, where q is the number of equality constraints, and q <= n. Equality constraints must be linearly independent with rank(Aeq) = q.

If your problem has no equality constraints, use [].

beq — Right-hand side of equality constraints

column vector of length q

Right-hand side of equality constraints, specified as a column vector of length q.

If your problem has no equality constraints, use zeros(0,1).

iA0 — Initial active inequalities

logical vector of length m

Initial active inequalities, where the equal portion of the inequality is true, specified as a logical vector of length m according to the following:

• If your problem has no inequality constraints, use false(0,1).

• For a cold start, false(m,1).

• For a warm start, set iA0(i) == true to start the algorithm with the ith inequality constraint active. Use the optional output argument iA from a previous solution to specify iA0 in this way. If both iA0(i) and iA0(j) are true, then rows i and j of A should be linearly independent. Otherwise, the solution can fail with status = -2.

options — Option set for mpcqpsolver

structure

Option set for mpcqpsolver, specified as a structure created using mpcqpsolverOptions.

Output Arguments

x — Optimal solution to the QP problem

column vector

Optimal solution to the QP problem, returned as a column vector of length n. mpcqpsolver always returns a value for x. To determine whether the solution is optimal or feasible, check the solution status.

status — Solution validity indicator

positive integer | 0 | -1 | -2

Solution validity indicator, returned as an integer according to the following:

Value	Description
> 0	x is optimal. status represents the number of iterations performed during optimization.
0	The maximum number of iterations was reached. The solution, x, may be suboptimal or infeasible.
-1	The problem appears to be infeasible, that is, the constraint $$Ax\ge b$$ cannot be satisfied.
-2	An unrecoverable numerical error occurred.

iA — Active inequalities

logical vector of length m

Active inequalities, where the equal portion of the inequality is true, returned as a logical vector of length m. If iA(i) == true, then the ith inequality is active for the solution x.

Use iA to warm start a subsequent mpcqpsolver solution.

lambda — Lagrange multipliers

structure

Lagrange multipliers, returned as a structure with the following fields:

Field	Description
ineqlin	Multipliers of the inequality constraints, returned as a vector of length n. When the solution is optimal, the elements of ineqlin are nonnegative.
eqlin	Multipliers of the equality constraints, returned as a vector of length q. There are no sign restrictions in the optimal solution.

Tips

• The KWIK algorithm requires that the Hessian matrix, H, be positive definite. When calculating Linv, use:

  [L, p] = chol(H,'lower');

  If p = 0, then H is positive definite. Otherwise, p is a positive integer.

• mpcqpsolver provides access to the QP solver used by Model Predictive Control Toolbox™ software. Use this command to solve QP problems in your own custom MPC applications.

Algorithms

mpcqpsolver solves the QP problem using an active-set method, the KWIK algorithm, based on [1]. For more information, see QP Solvers.

References

[1] Schmid, C., and L.T. Biegler. ‘Quadratic Programming Methods for Reduced Hessian SQP’. Computers & Chemical Engineering 18, no. 9 (September 1994): 817–32. https://doi.org/10.1016/0098-1354(94)E0001-4.

Extended Capabilities

C/C++ Code Generation

Generate C and C++ code using MATLAB® Coder™.

Usage notes and limitations:

• You can use mpcqpsolver as a general-purpose QP solver that supports code generation. Create the function myCode that uses mpcqpsolver.

  function [out1,out2] = myCode(in1,in2)
  %#codegen
  ...
  [x,status] = mpcqpsolver(Linv,f,A,b,Aeq,Beq,iA0,options);
  ...

  Generate C code with MATLAB^® Coder™.

  func = 'myCode';
  cfg = coder.config('mex'); % or 'lib', 'dll'
  codegen('-config',cfg,func,'-o',func);

• For code generation, use the same precision for all real inputs, including options. Configure the precision as 'double' or 'single' using mpcqpsolverOptions.

Version History

Introduced in R2015b

R2020a: mpcqpsolver will be removed

mpcqpsolver will be removed in a future release. Use mpcActiveSetSolver instead. There are differences between these functions that require updates to your code.

Update Code

The following differences require updates to your code:

• For mpcActiveSetSolver, you define inequality constraints in the form Ax≤b. Previously, for mpcqpsolver, you defined inequality constraints in the form Ax≥b

• For mpcActiveSetSolver, you specify solver options with a structure created using the mpcActiveSetOptions function. Previously, for mpcqpsolver, you created an option structure using the mpcqpsolverOptions function. These option structures contain the same options, though some option names have changed.

• By default, you pass the Hessian matrix to mpcActiveSetSolver. Previously, you passed the inverse of lower-triangular Cholesky decomposition (Linv) of the Hessian matrix to mpcqpsolver. To continue to use Linv, set the UseHessianAsInput field of the structure returned by mpcActiveSetSolver to false.

• When your QP problem has either no inequality constraints or no equality constraints, the corresponding A or Aeq input argument to mpcActiveSetSolver must be zeros(0,n), where n is the number of decision variables. Previously, for mpcqpsolver, you specified these input arguments as [].

This table shows some typical usages of mpcqpsolver and how to update your code to use mpcActiveSetSolver instead.

Not Recommended	Recommended
opt = mpcqpsolverOptions; [x,status] = mpcqpsolver(Linv,f,A,b,... Aeq,beq,iA0,opt);	opt = mpcActiveSetOptions; opt.UseHessianAsInput = false; [x,status] = mpcActiveSetSolver(Linv,f,... -A,-b,Aeq,beq,iA0,opt); Alternatively, you can use the Hessian matrix, H. opt = mpcActiveSetOptions; [x,status] = mpcActiveSetSolver(H,f,... -A,-b,Aeq,beq,iA0,opt);
opt = mpcqpsolverOptions('single'); [x,status] = mpcqpsolver(Linv,f,A,b,... Aeq,beq,iA0,opt);	opt = mpcActiveSetOptions('single'); opt.UseHessianAsInput = false; [x,status] = mpcActiveSetSolver(Linv,f,... -A,-b,Aeq,beq,iA0,opt);
opt = mpcqpsolverOptions; opt.MaxIter = 300; opt.FeasibilityTol = 1e-5; [x,status] = mpcqpsolver(Linv,f,A,b,... Aeq,beq,iA0,opt);	opt = mpcActiveSetOptions; opt.UseHessianAsInput = false; opt.MaxIterations = 300; opt.ContraintTolerance = 1e-5; [x,status] = mpcActiveSetSolver(Linv,f,... -A,-b,Aeq,beq,iA0,opt);
[x,status] = mpcqpsolver(Linv,f,[],... zeros(0,1),Aeq,beq,iA0,opt);	n = length(f); opt.UseHessianAsInput = false; [x,status] = mpcActiveSetSolver(Linv,f,... zeros(0,n),zeros(0,1),Aeq,beq,iA0,opt);
[x,status] = mpcqpsolver(Linv,f,A,b,... [],zeros(0,1),iA0,opt);	n = length(f); opt.UseHessianAsInput = false; [x,status] = mpcActiveSetSolver(Linv,f,... -A,-b,zeros(0,n),zeros(0,1),iA0,opt);