具有二次约束的线性或二次目标

此示例说明如何求解具有线性或二次目标和二次不等式约束的优化问题。该示例生成并使用目标函数和约束函数的梯度和黑塞函数。

二次约束问题
目标函数
约束函数
数值示例
黑塞矩阵
解
提供黑塞函数时的效率
二次等式约束的扩展

二次约束问题

假设您的问题具有以下形式

$\min_{x} (\frac{1}{2} x^{T} Q x + f^{T} x + c)$

约束条件为

$\frac{1}{2} x^{T} H_{i} x + k_{i}^{T} x + d_{i} \leq 0,$

其中 1 ≤ i ≤ m。假设有至少一个 H_i 为非零值；否则，您可以使用 quadprog 或 linprog 来求解此问题。使用非零 H_i 时，约束是非线性的，这意味着根据优化决策表，fmincon 是合适的求解器。

此示例假设二次矩阵是对称的，不失其普遍性。您可以将非对称 H（或 Q）矩阵替换为等效的对称版本 (H + H^T)/2。

如果 x 包含 N 个分量，则 Q 和 H_i 是 N×N 矩阵，f 和 k_i 是 N×1 向量，c 和 d_i 是标量。

目标函数

使用 fmincon 语法表示问题。假设 x 和 f 是列向量。（当初始向量 x0 是列向量时，x 是列向量。）

function [y,grady] = quadobj(x,Q,f,c)
y = 1/2*x'*Q*x + f'*x + c;
if nargout > 1
    grady = Q*x + f;
end

约束函数

为了一致性和易于索引，将每个二次约束矩阵放在一个元胞数组中。同样，将线性项和常数项也放在元胞数组中。

function [y,yeq,grady,gradyeq] = quadconstr(x,H,k,d)
jj = length(H); % jj is the number of inequality constraints
y = zeros(1,jj);
for i = 1:jj
    y(i) = 1/2*x'*H{i}*x + k{i}'*x + d{i};
end
yeq = [];
    
if nargout > 2
    grady = zeros(length(x),jj);
    for i = 1:jj
        grady(:,i) = H{i}*x + k{i};
    end
end
gradyeq = [];

数值示例

假设您有以下问题。

Q = [3,2,1;
     2,4,0;
     1,0,5];
f = [-24;-48;-130];
c = -2;

rng default % For reproducibility
% Two sets of random quadratic constraints:
H{1} = gallery('randcorr',3); % Random positive definite matrix
H{2} = gallery('randcorr',3);
k{1} = randn(3,1);
k{2} = randn(3,1);
d{1} = randn;
d{2} = randn;

黑塞矩阵

创建一个黑塞函数。拉格朗日函数的黑塞矩阵由以下方程给出

$\nabla_{x x}^{2} L (x, λ) = \nabla^{2} f (x) + \sum λ_{i} \nabla^{2} c_{i} (x) + \sum λ_{i} \nabla^{2} c e q_{i} (x) .$

fmincon 计算一组近似拉格朗日乘数 λ_i，并将它们放在一个结构体中。要包含黑塞函数，请使用以下函数。

function hess = quadhess(x,lambda,Q,H)
hess = Q;
jj = length(H); % jj is the number of inequality constraints
for i = 1:jj
    hess = hess + lambda.ineqnonlin(i)*H{i};
end

解

使用 fmincon interior-point 算法可最高效地求解问题。该算法接受您提供的黑塞函数。设置以下选项。

options = optimoptions(@fmincon,'Algorithm','interior-point',...
    'SpecifyObjectiveGradient',true,'SpecifyConstraintGradient',true,...
    'HessianFcn',@(x,lambda)quadhess(x,lambda,Q,H));

调用 fmincon 以求解问题。

fun = @(x)quadobj(x,Q,f,c);
nonlconstr = @(x)quadconstr(x,H,k,d);
x0 = [0;0;0]; % Column vector
[x,fval,eflag,output,lambda] = fmincon(fun,x0,...
    [],[],[],[],[],[],nonlconstr,options);

检查拉格朗日乘数。

lambda.ineqnonlin

ans =

   12.8412
   39.2337

两个非线性不等式乘数均为非零，因此两个二次约束在解处都为活动约束。

提供黑塞函数时的效率

使用梯度和黑塞函数的内点算法很高效。查看函数计算的次数。

output

output = 

         iterations: 9
          funcCount: 10
    constrviolation: 0
           stepsize: 5.3547e-04
          algorithm: 'interior-point'
      firstorderopt: 1.5851e-05
       cgiterations: 0
            message: 'Local minimum found that satisfies the constraints.

Optimization compl...'

fmincon 只需 10 次函数计算即可求解问题。

将此结果与不使用黑塞函数时的解进行比较。

options.HessianFcn = [];
[x2,fval2,eflag2,output2,lambda2] = fmincon(fun,[0;0;0],...
    [],[],[],[],[],[],nonlconstr,options);
output2

output2 = 

         iterations: 17
          funcCount: 22
    constrviolation: 0
           stepsize: 2.8475e-04
          algorithm: 'interior-point'
      firstorderopt: 1.7680e-05
       cgiterations: 0
            message: 'Local minimum found that satisfies the constraints.

Optimization compl...'

在这种情况下，fmincon 需要大约两倍的迭代次数和函数计算次数。解都同样在容差范围内。