Esta página aún no se ha traducido para esta versión. Puede ver la versión más reciente de esta página en inglés.

Calcula gradientes y hessianosSymbolic Math Toolbox

Si tiene una licencia, puede calcular fácilmente gradientes analíticos y hessianos para funciones objetivas y de restricción.Symbolic Math Toolbox™ Existen dos funciones relevantes:Symbolic Math Toolbox

  • genera el degradado de una función escalar y genera una matriz de los derivados parciales de una función vectorial.jacobian Así, por ejemplo, se puede obtener la matriz Hessiana, los segundos derivados de la función objetiva, aplicando al degradado.jacobian Parte de este ejemplo muestra cómo utilizar para generar degradados simbólicos y hessianos de funciones objetivas y de restricción.jacobian

  • genera una función anónima o un archivo que calcula los valores de una expresión simbólica.matlabFunction En este ejemplo se muestra cómo usar para generar archivos que evalúan la función de objetivo y restricción y sus derivados en puntos arbitrarios.matlabFunction

Consideremos el problema de la electroestática de colocar 10 electrones en un cuerpo conductor. Los electrones se arreglaran para minimizar su energía potencial total, sujeto a la restricción de acostarse dentro del cuerpo. Es bien sabido que todos los electrones estarán en el límite del cuerpo como mínimo. Los electrones son indistinguibles, por lo que no hay un mínimo único para este problema (permuting los electrones en una solución da otra solución válida). Este ejemplo fue inspirado por Dolan, Moré y Munson.[58]

Este ejemplo es un organismo conductor definido por las siguientes desigualdades:

z|x||y|(1)
x2+y2+(z+1)21.(2)

Este cuerpo se ve como una pirámide en una esfera.

 Código para generar la figura

Hay un pequeño hueco entre las superficies superior e inferior de la figura. Este es un artefacto de la rutina de trazado general utilizada para crear la figura. Esta rutina borra cualquier parche rectangular en una superficie que toque la otra superficie.

La sintaxis y las estructuras de los dos conjuntos de funciones de Toolbox difieren. En particular, las variables simbólicas son escalares reales o complejos, pero las funciones pasan argumentos vectoriales.Optimization Toolbox™ Así que hay varios pasos a seguir para generar simbólicamente la función objetiva, las restricciones, y todos sus derivados requeridos, en una forma adecuada para el algoritmo de punto interior de:fmincon

Para ver la eficiencia en el uso de gradientes y hessianos, ver.Comparar con optimización sin gradientes y hessianos Para obtener un enfoque basado en problemas de este problema sin utilizar información derivada, consulte.Optimización no lineal electrostática restringida, basada en problemas

Cree las variables

Genere un vector simbólico como un vector de 30 por 1 compuesto de variables simbólicas reales, entre 1 y 10, y entre 1 y 3.xxijij Estas variables representan las tres coordenadas del electrón: corresponde a la coordenada, corresponde a la coordenada, y corresponde a la coordenada.ixi1xxi2yxi3z

x = cell(3, 10); for i = 1:10     for j = 1:3         x{j,i} = sprintf('x%d%d',i,j);     end end x = x(:); % now x is a 30-by-1 vector x = sym(x, 'real');

El vector es:x

x x =    x11   x12   x13   x21   x22   x23   x31   x32   x33   x41   x42   x43   x51   x52   x53   x61   x62   x63   x71   x72   x73   x81   x82   x83   x91   x92   x93  x101  x102  x103

Incluya las restricciones lineales

Escriba las restricciones lineales en,Ecuación 1

≤ – | | – | |,zxy

como un conjunto de cuatro desigualdades lineales para cada electrón:

– – ≤ 0 – + ≤ 0 + – ≤ 0 + + ≤ 0xi3xi1xi2
xi3xi1xi2
xi3xi1xi2
xi3xi1xi2

Por lo tanto, hay un total de 40 desigualdades lineales para este problema.

Escribir las desigualdades de una manera estructurada:

B = [1,1,1;-1,1,1;1,-1,1;-1,-1,1];  A = zeros(40,30); for i=1:10     A(4*i-3:4*i,3*i-2:3*i) = B; end  b = zeros(40,1);

Se puede ver que representa las desigualdades:A*x ≤ b

A*x   ans =      x11 + x12 + x13     x12 - x11 + x13     x11 - x12 + x13     x13 - x12 - x11     x21 + x22 + x23     x22 - x21 + x23     x21 - x22 + x23     x23 - x22 - x21     x31 + x32 + x33     x32 - x31 + x33     x31 - x32 + x33     x33 - x32 - x31     x41 + x42 + x43     x42 - x41 + x43     x41 - x42 + x43     x43 - x42 - x41     x51 + x52 + x53     x52 - x51 + x53     x51 - x52 + x53     x53 - x52 - x51     x61 + x62 + x63     x62 - x61 + x63     x61 - x62 + x63     x63 - x62 - x61     x71 + x72 + x73     x72 - x71 + x73     x71 - x72 + x73     x73 - x72 - x71     x81 + x82 + x83     x82 - x81 + x83     x81 - x82 + x83     x83 - x82 - x81     x91 + x92 + x93     x92 - x91 + x93     x91 - x92 + x93     x93 - x92 - x91  x101 + x102 + x103  x102 - x101 + x103  x101 - x102 + x103  x103 - x102 - x101

Cree las restricciones no lineales, sus degradados y los hessianos

Las restricciones no lineales en,Ecuación 2

x2+y2+(z+1)21,

también están estructurados. Genere las restricciones, sus gradientes y los hessianos de la siguiente manera:

c = sym(zeros(1,10)); i = 1:10; c = (x(3*i-2).^2 + x(3*i-1).^2 + (x(3*i)+1).^2 - 1).';  gradc = jacobian(c,x).'; % .' performs transpose  hessc = cell(1, 10); for i = 1:10     hessc{i} = jacobian(gradc(:,i),x); end

El vector de restricción es un vector de fila, y el degradado de se representa en la columna TH de la matriz.cc(i)igradc Esta es la forma correcta, como se describe en.Restricciones no lineales

Las matrices de hessian,..., son cuadradas y simétricas.hessc{1}hessc{10} Es mejor almacenarlos en una matriz de celdas, como se hace aquí, que en variables separadas como.hessc1, ..., hesssc10

Utilice la sintaxis para transponer..' La sintaxis significa transposición conjugada, que tiene diferentes derivados simbólicos.'

Crear la función objetivo, su gradiente y hessian

La función objetiva, energía potencial, es la suma de las inversas de las distancias entre cada par de electrones:

energy=i<j1|xixj|.

La distancia es la raíz cuadrada de la suma de los cuadrados de las diferencias en los componentes de los vectores.

Calcule la energía, su gradiente y su hessian de la siguiente manera:

energy = sym(0); for i = 1:3:25     for j = i+3:3:28         dist = x(i:i+2) - x(j:j+2);         energy = energy + 1/sqrt(dist.'*dist);     end end  gradenergy = jacobian(energy,x).';  hessenergy = jacobian(gradenergy,x);

Cree el archivo de función objetivo

La función objetiva debe tener dos salidas, y.energygradenergy Coloque ambas funciones en un vector al llamar para reducir el número de subexpresiones que genera y devolver el degradado sólo cuando la función de llamada (en este caso) solicita ambas salidas.matlabFunctionmatlabFunctionfmincon Este ejemplo muestra la colocación de los archivos resultantes en la carpeta actual. Por supuesto, puede colocarlos en cualquier lugar que quiera, siempre y cuando la carpeta esté en la ruta de MATLAB.

currdir = [pwd filesep]; % You may need to use currdir = pwd  filename = [currdir,'demoenergy.m']; matlabFunction(energy,gradenergy,'file',filename,'vars',{x});

Esta sintaxis provoca que se devuelva como la primera salida y como la segunda.matlabFunctionenergygradenergy También toma un vector de entrada único en lugar de una lista de entradas,...,.{x}x11x103

El archivo resultante contiene, en parte, las siguientes líneas u otras similares:demoenergy.m

function [energy,gradenergy] = demoenergy(in1) %DEMOENERGY %   [ENERGY,GRADENERGY] = DEMOENERGY(IN1) ... x101 = in1(28,:); ... energy = 1./t140.^(1./2) + ...; if nargout > 1     ...     gradenergy = [(t174.*(t185 - 2.*x11))./2 - ...]; end

Esta función tiene la forma correcta para una función objetiva con un degradado; Ver.Escribir funciones de objetivo escalar

Cree el archivo de función de restricción

Genere la función de restricción no lineal y colóque en el formato correcto.

filename = [currdir,'democonstr.m']; matlabFunction(c,[],gradc,[],'file',filename,'vars',{x},...     'outputs',{'c','ceq','gradc','gradceq'});

El archivo resultante contiene, en parte, las siguientes líneas u otras similares:democonstr.m

function [c,ceq,gradc,gradceq] = democonstr(in1) %DEMOCONSTR %    [C,CEQ,GRADC,GRADCEQ] = DEMOCONSTR(IN1) ... x101 = in1(28,:); ... c = [t417.^2 + ...]; if nargout > 1     ceq = []; end if nargout > 2     gradc = [2.*x11,...]; end if nargout > 3     gradceq = []; end 

Esta función tiene la forma correcta para una función de restricción con un degradado; Ver.Restricciones no lineales

Genere los archivos de hessian

Para generar el hessian del Lagrangio para el problema, primero generar archivos para la energía hessian y para la restricción de los hessianos.

El Hessiano de la función objetiva,, es una expresión simbólica muy grande, que contiene más de 150.000 símbolos, como la evaluación de espectáculos.hessenergysize(char(hessenergy)) Así que se necesita una cantidad considerable de tiempo para correr.matlabFunction(hessenergy)

Para generar un archivo, ejecute las dos líneas siguientes:hessenergy.m

filename = [currdir,'hessenergy.m']; matlabFunction(hessenergy,'file',filename,'vars',{x});

Por el contrario, los hessianos de las funciones de restricción son pequeños y rápidos de calcular:

for i = 1:10     ii = num2str(i);     thename = ['hessc',ii];     filename = [currdir,thename,'.m'];     matlabFunction(hessc{i},'file',filename,'vars',{x}); end

Después de generar todos los archivos para el objetivo y las restricciones, Póntelos junto con los multiplicadores de Lagrange apropiados en un archivo de la siguiente manera:hessfinal.m

function H = hessfinal(X,lambda) % % Call the function hessenergy to start H = hessenergy(X);  % Add the Lagrange multipliers * the constraint Hessians H = H + hessc1(X) * lambda.ineqnonlin(1); H = H + hessc2(X) * lambda.ineqnonlin(2); H = H + hessc3(X) * lambda.ineqnonlin(3); H = H + hessc4(X) * lambda.ineqnonlin(4); H = H + hessc5(X) * lambda.ineqnonlin(5); H = H + hessc6(X) * lambda.ineqnonlin(6); H = H + hessc7(X) * lambda.ineqnonlin(7); H = H + hessc8(X) * lambda.ineqnonlin(8); H = H + hessc9(X) * lambda.ineqnonlin(9); H = H + hessc10(X) * lambda.ineqnonlin(10);

Ejecute la optimización

Iniciar la optimización con los electrones distribuidos aleatoriamente en una esfera de radio 1/2 centrado en [0, 0, – 1]:

rng default % for reproducibility Xinitial = randn(3,10); % columns are normal 3-D vectors for j=1:10     Xinitial(:,j) = Xinitial(:,j)/norm(Xinitial(:,j))/2;     % this normalizes to a 1/2-sphere end Xinitial(3,:) = Xinitial(3,:) - 1; % center at [0,0,-1] Xinitial = Xinitial(:); % Convert to a column vector

Establezca las opciones para utilizar el algoritmo de punto interior, y para utilizar gradientes y el hessian:

options = optimoptions(@fmincon,'Algorithm','interior-point','SpecifyObjectiveGradient',true,...     'SpecifyConstraintGradient',true,'HessianFcn',@hessfinal,'Display','final');

Llamar:fmincon

[xfinal fval exitflag output] = fmincon(@demoenergy,Xinitial,...     A,b,[],[],[],[],@democonstr,options);

La solución toma 19 iteraciones y solo 28 evaluaciones de función:

xfinal,fval,exitflag,output.iterations,output.funcCount
xfinal =     -0.0317     0.0317    -1.9990     0.6356    -0.6356    -1.4381     0.0000    -0.0000    -0.0000     0.0000    -1.0000    -1.0000     1.0000    -0.0000    -1.0000    -1.0000    -0.0000    -1.0000     0.6689     0.6644    -1.3333    -0.6667     0.6667    -1.3333     0.0000     1.0000    -1.0000    -0.6644    -0.6689    -1.3333   fval =     34.1365   exitflag =       1   ans =      19   ans =      28

A pesar de que las posiciones iniciales de los electrones eran aleatorias, las posiciones finales son casi simétricas:

 Código para generar la figura

Comparar con optimización sin gradientes y hessianos

El uso de gradientes y hessianos hace que la optimización se ejecute más rápido y con mayor precisión. Para comparar con la misma optimización sin utilizar información de degradado o de hessian, establezca las opciones para no utilizar gradientes y hessianos:

options = optimoptions(@fmincon,'Algorithm','interior-point',...     'Display','final');
[xfinal2 fval2 exitflag2 output2] = fmincon(@demoenergy,Xinitial,...     A,b,[],[],[],[],@democonstr,options);

La salida muestra que encontró un mínimo equivalente, pero tomó más iteraciones y muchas más evaluaciones de funciones para hacerlo.fmincon

xfinal2,fval2,exitflag2,output2.iterations,output2.funcCount
xfinal2 =      0.0000     1.0000    -1.0000     0.6689    -0.6644    -1.3334    -0.6644     0.6689    -1.3334     0.0000    -1.0000    -1.0000     0.6357     0.6357    -1.4380    -0.0317    -0.0317    -1.9990     1.0000     0.0000    -1.0000    -1.0000     0.0000    -1.0000     0.0000     0.0000    -0.0000    -0.6667    -0.6667    -1.3334   fval2 =     34.1365   exitflag2 =       1   ans =      77   ans =          2435

En esta ejecución, el número de evaluaciones de función (in) es 2435 comparado con 28 (in) cuando se utilizan gradientes y hessian.output2.funcCountoutput.funcCount

Desactive las hipótesis de variables simbólicas

Las variables simbólicas en este ejemplo tienen la suposición, en el espacio de trabajo del motor simbólico, de que son reales. Para borrar esta suposición del espacio de trabajo del motor simbólico, no es suficiente eliminar las variables. Desactive las suposiciones de variables mediante:syms

syms x

Compruebe que las suposiciones están vacías.

assumptions(x)
ans =   Empty sym: 1-by-0

Temas relacionados