Main Content

Optimización no lineal electrostática restringida, basada en problemas

Abrir script en vivo

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

Para obtener un ejemplo basado en un solucionador equivalente mediante el cuadro de herramientas de matemática simbólica™, consulte.Calcula gradientes y hessianosSymbolic Math Toolbox

Problema Geometry

Este ejemplo implica un organismo conductor definido por las siguientes desigualdades. Para cada electrón con coordenadas

<math display="block">
<mrow>
<mo stretchy="false">(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
<mo>,</mo>
<mi>z</mi>
<mo stretchy="false">)</mo>
</mrow>
</math>
,

<math display="block">
<mrow>
<mtable columnalign="left">
<mtr>
<mtd>
<mrow>
<mi>z</mi>
<mo>≤</mo>
<mo>-</mo>
<mo>|</mo>
<mi>x</mi>
<mo>|</mo>
<mo>-</mo>
<mo>|</mo>
<mi>y</mi>
<mo>|</mo>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<msup>
<mrow>
<mi>x</mi>
</mrow>
<mrow>
<mn>2</mn>
</mrow>
</msup>
<mo>+</mo>
<msup>
<mrow>
<mi>y</mi>
</mrow>
<mrow>
<mn>2</mn>
</mrow>
</msup>
<mo>+</mo>
<mo stretchy="false">(</mo>
<mi>z</mi>
<mo>+</mo>
<mn>1</mn>
<msup>
<mrow>
<mo stretchy="false">)</mo>
</mrow>
<mrow>
<mn>2</mn>
</mrow>
</msup>
<mo>≤</mo>
<mn>1</mn>
<mo>.</mo>
</mrow>
</mtd>
</mtr>
</mtable>
</mrow>
</math>

Estas restricciones forman un cuerpo que se ve como una pirámide en una esfera. Para ver el cuerpo, escriba el código siguiente.

[X,Y] = meshgrid(-1:.01:1); Z1 = -abs(X) - abs(Y); Z2 = -1 - sqrt(1 - X.^2 - Y.^2); Z2 = real(Z2); W1 = Z1; W2 = Z2; W1(Z1 < Z2) = nan; % only plot points where Z1 > Z2 W2(Z1 < Z2) = nan; % only plot points where Z1 > Z2 hand = figure; % handle to the figure, since we'll plot more later set(gcf,'Color','w') % white background surf(X,Y,W1,'LineStyle','none'); hold on surf(X,Y,W2,'LineStyle','none'); view(-44,18)

Existe una ligera brecha entre las superficies superior e inferior de la figura. Esta brecha es un artefacto de la rutina de trazado general utilizada para crear la figura. La rutina borra cualquier parche rectangular en una superficie que toque la otra superficie.

Definir variables de problema

El problema tiene diez electrones. Las restricciones dan límites a cada

<math display="block">
<mrow>
<mi>x</mi>
</mrow>
</math>
Y
<math display="block">
<mrow>
<mi>y</mi>
</mrow>
</math>
 valor de – 1 a 1, y el
<math display="block">
<mrow>
<mi>z</mi>
</mrow>
</math>
 valor de – 2 a 0. Defina las variables para el problema.

N = 10; x = optimvar('x',N,'LowerBound',-1,'UpperBound',1); y = optimvar('y',N,'LowerBound',-1,'UpperBound',1); z = optimvar('z',N,'LowerBound',-2,'UpperBound',0); elecprob = optimproblem;

Definir restricciones

El problema tiene dos tipos de restricciones. La primera, una restricción esférica, es una simple desigualdad polinómica para cada electrón por separado. Defina esta restricción esférica.

elecprob.Constraints.spherec = (x.^2 + y.^2 + (z+1).^2) <= 1;

El comando de restricción anterior crea un vector de diez restricciones. Ver el vector de restricción utilizando.showconstr

showconstr(elecprob.Constraints.spherec)
  ((x.^2 + y.^2) + (z + 1).^2) <= arg_RHS    where:      arg2 = 1;     arg1 = arg2([1 1 1 1 1 1 1 1 1 1]);     arg_RHS = arg1(:);

El segundo tipo de restricción en el problema es lineal. Puede expresar las restricciones lineales de diferentes maneras. Por ejemplo, puede utilizar la función para representar una restricción de valor absoluto.abs Para expresar las restricciones de esta manera, escriba una función de MATLAB y conviértalo en una expresión utilizando.fcn2optimexpr Para un enfoque diferente, escriba la restricción de valor absoluto como cuatro desigualdades lineales. Cada comando de restricción devuelve un vector de diez restricciones.

elecprob.Constraints.plane1 = z <= -x-y; elecprob.Constraints.plane2 = z <= -x+y; elecprob.Constraints.plane3 = z <= x-y; elecprob.Constraints.plane4 = z <= x+y;

Definir función objetiva

La función objetiva es la energía potencial del sistema, que es una suma sobre cada par de electrones de la inversa de sus distancias:

<math display="block">
<mrow>
<mrow>
<mstyle mathvariant="normal">
<mrow>
<mi>e</mi>
<mi>n</mi>
<mi>e</mi>
<mi>r</mi>
<mi>g</mi>
<mi>y</mi>
</mrow>
</mstyle>
</mrow>
<mo>=</mo>
<munder>
<mrow>
<mo>∑</mo>
</mrow>
<mrow>
<mi>i</mi>
<mo><</mo>
<mi>j</mi>
</mrow>
</munder>
<mfrac>
<mrow>
<mn>1</mn>
</mrow>
<mrow>
<mo stretchy="false">‖</mo>
<mrow>
<mstyle mathvariant="normal">
<mrow>
<mi>e</mi>
<mi>l</mi>
<mi>e</mi>
<mi>c</mi>
<mi>t</mi>
<mi>r</mi>
<mi>o</mi>
<mi>n</mi>
</mrow>
</mstyle>
</mrow>
<mo stretchy="false">(</mo>
<mi>i</mi>
<mo stretchy="false">)</mo>
<mo>-</mo>
<mrow>
<mstyle mathvariant="normal">
<mrow>
<mi>e</mi>
<mi>l</mi>
<mi>e</mi>
<mi>c</mi>
<mi>t</mi>
<mi>r</mi>
<mi>o</mi>
<mi>n</mi>
</mrow>
</mstyle>
</mrow>
<mo stretchy="false">(</mo>
<mi>j</mi>
<mo stretchy="false">)</mo>
<mo stretchy="false">‖</mo>
</mrow>
</mfrac>
<mo>.</mo>
</mrow>
</math>

Defina la función objetiva como una expresión de optimización.

energy = optimexpr(1); for ii = 1:(N-1)     for jj = (ii+1):N         tempe = (x(ii) - x(jj))^2 + (y(ii) - y(jj))^2 + (z(ii) - z(jj))^2;         energy = energy + tempe^(-1/2);     end end elecprob.Objective = energy;

Ejecute Optimization

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 x0 = randn(N,3); for ii=1:N     x0(ii,:) = x0(ii,:)/norm(x0(ii,:))/2;     x0(ii,3) = x0(ii,3) - 1; end init.x = x0(:,1); init.y = x0(:,2); init.z = x0(:,3);

Resuelva el problema llamando.solve

[sol,fval,exitflag,output] = solve(elecprob,init)
Local minimum found that satisfies the constraints.  Optimization completed because the objective function is non-decreasing in  feasible directions, to within the value of the optimality tolerance, and constraints are satisfied to within the value of the constraint tolerance. 
sol = struct with fields:
    x: [10x1 double]
    y: [10x1 double]
    z: [10x1 double]


fval = 34.1365 
exitflag =      OptimalSolution  
output = struct with fields:
         iterations: 95
          funcCount: 2990
    constrviolation: 0
           stepsize: 6.4513e-07
          algorithm: 'interior-point'
      firstorderopt: 1.2800e-06
       cgiterations: 0
            message: '...'
             solver: 'fmincon'

Ver solución

Trace la solución como puntos en el cuerpo conductor.

figure(hand) plot3(sol.x,sol.y,sol.z,'r.','MarkerSize',25) hold off

Los electrones se distribuyen de manera bastante uniforme en el contorno de la restricción. Los electrones se encuentran principalmente en los bordes y el punto piramidal.

Referencia

[1] Dolan, Elizabeth D., Jorge J. Moré, y Todd S. Munson. "Benchmarking software de optimización con COPS 3,0." Informe técnico del laboratorio nacional de Argonne ANL/MCS-TM-273, febrero de 2004.

Temas relacionados

Documentación de Optimization Toolbox

Soporte

8 MATLAB Cheat Sheets for Data Science

8 MATLAB Cheat Sheets for Data Science

Download now