Objetivo lineal o cuadrático con restricciones cuadráticas
Este ejemplo muestra cómo resolver un problema de optimización que tiene un objetivo lineal o cuadrático y restricciones de desigualdad cuadráticas. El ejemplo genera y usa el gradiente y la matriz hessiana de la función objetivo y de la función de restricción.
Problema cuadrático restringido
Suponga que su problema tiene el formato
sujeto a
donde 1 ≤ i ≤ m. Asuma que al menos una Hi es distinta de cero; de lo contrario, puede usar quadprog
o linprog
para resolver este problema. Con Hi distintas de cero, las restricciones son no lineales, lo que significa que fmincon
es el solver adecuado de acuerdo con Tabla de decisiones de optimización.
El ejemplo da por supuesto que las matrices cuadráticas son simétricas sin pérdida de generalidad. Puede reemplazar una matriz H (o Q) no simétrica con una versión equivalente simetrizada (H + HT)/2.
Si x tiene N componentes, Q y Hi son matrices de N por N, f y ki son vectores de N por 1 y c y di son escalares.
Función objetivo
Formule el problema utilizando la sintaxis de fmincon
. Asuma que x
y f
son vectores columna. (x
es un vector columna cuando el vector inicial x0
es un vector columna).
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
Función de restricción
Por motivos de consistencia y para una indexación fácil, coloque cada matriz de restricción cuadrática en un arreglo de celdas. De forma similar, coloque el término lineal y el término constante en arreglos de celdas.
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 = [];
Ejemplo numérico
Suponga que tiene el siguiente problema.
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;
Matriz hessiana
Cree una función de matriz hessiana. La ecuación proporciona la matriz hessiana del lagrangiano
fmincon
calcula un conjunto aproximado de multiplicadores de Lagrange λi y los agrupa en una estructura. Para incluir la matriz hessiana, utilice la siguiente función.
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
Solución
Utilice el algoritmo interior-point
de fmincon
para resolver el problema con la máxima eficiencia. Este algoritmo acepta que proporcione una función de matriz hessiana. Configure estas opciones.
options = optimoptions(@fmincon,'Algorithm','interior-point',... 'SpecifyObjectiveGradient',true,'SpecifyConstraintGradient',true,... 'HessianFcn',@(x,lambda)quadhess(x,lambda,Q,H));
Llame a fmincon
para resolver el problema.
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);
Examine los multiplicadores de Lagrange.
lambda.ineqnonlin
ans = 12.8412 39.2337
Ambos multiplicadores de desigualdad no lineal son distintos de cero, por lo que ambas restricciones cuadráticas están activas en la solución.
Eficiencia cuando se proporciona una matriz hessiana
El algoritmo interior-point con gradientes y una matriz hessiana es eficiente. Visualice el número de evaluaciones de función.
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
necesita solo 10 evaluaciones de función para resolver el problema.
Compare este resultado con la solución sin la matriz hessiana.
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...'
En este caso, fmincon
necesita casi el doble de iteraciones y evaluaciones de función. Las soluciones son las mismas que dentro de las tolerancias.
Extensión para las restricciones de igualdad cuadráticas
Si también tiene restricciones de igualdad cuadráticas, puede usar básicamente la misma técnica. El problema es el mismo, con las restricciones adicionales
Reformule las restricciones para usar las variables Ji, pi y qi. La estructura de lambda.eqnonlin(i)
tiene los multiplicadores de Lagrange para restricciones de igualdad.