Cómo calcular el mínimo de una función con restricciones

La función "fmincon" se utiliza para resolver problemas de este tipo:

calcular x, minimizando fun(x)

sujeto a las siguientes restricciones:

• g(x) <= 0
• geq(x) = 0
• A·x - B <= 0
• C·x - D = 0
• xmin <= x <= xmax

Dicha función se utiliza de las siguientes maneras:

• fmincon(@(x) fun(x), x0, A,B)
• fmincon(@(x) fun(x), x0, A,B,C,D,xmin,xmax)
• fmincon(@(x) fun(x), x0, A,B,C,D,xmin,xmax,@(x) funrest(x))

donde:

• x0 es el valor inicial a partir del que Matlab empieza a tantear.
• Si uno de los parámetros no existe, se pone [].
• fun(x) es la función a optimizar, y devuelve un escalar.
• las restricciones lineales se añaden con las matrices A, B, C y D.
• las restricciones no lineales se escriben dentro de una función ("funrest(x)"), que devuelve dos dos vectores, g(x) y geq(x). (El tamaño de los dos vectores corresponde con el número de restricciones de tipo "menor o igual que" e "igual que".)

---

Por ejemplo, queremos optimizar la función de Rosenbrock siguiente:

z = (1-x)² + 2·(y-x²)²

teniendo en cuenta las siguientes restricciones lineales:

• y + x -1 <= 0
• -y -1/3·x - 2/3 <= 0

y las siguientes restricciones no lineales:

• y + 0.94·x² - 0.14·x -3 <= 0
• -y +0.143·x² - 1.03 <= 0

---

Primero, creamos el fichero "f01.m" con la función a minimizar:

function z = f01(u)
x = u(1);
y = u(2);
z = (1-x).^2 + 2*(y-x.^2).^2;

A continuación, creamos el fichero "frest01.m" con las restricciones no lineales.

function [z zz] = fc01(x)

z(1) = x(2)+0.9375*x(1).^2-0.125*x(1)-3;
z(2) = -x(2)+0.1333*x(1).^2-1.0333;

zz = [];

En este caso, "z" son las restricciones de tipo "menor o igual que", mientras que "zz" son las restricciones de tipo "igual que". Como en este ejemplo, no hay restricciones del segundo tipo, la salida de "zz" es el conjunto vacío.

---

El mínimo de la función de Rosenbrock está en (1,1). En la siguiente figura se pueden ver las isolíneas de dicha función.

Para representar la anterior figura, lo hice con las siguientes líneas:

x = [-2:0.1:2];
y = [-1:0.1:3];
nx = length(x);
ny = length(y);
[X,Y] = meshgrid(x,y);
Z = zeros(nx,ny);
for i=1:1:nx
    for j=1:1:ny
        Z(i,j) = f01([X(i,j),Y(i,j)]);
    end
end

[u,f] = fminsearch(@(u) f01(u),[0 0]);

figure(1)
hold on
contour(X,Y,Z,[75])
plot(u(1),u(2),'+k')
set(gca,'XTick',[-2:1:2])
set(gca,'YTick',[-1:1:3])
set(gca,'Fontsize',10)
axis([-2 2 -1 3])
box on

---

A continuación, represento el área permitida a partir de las restricciones del ejemplo. Como se ve, el punto (1,1) está fuera de este área.

Para representar el área permitida, lo hice con las siguientes líneas:

t = linspace(-1.7476,-0.9785,100)';
r = -0.9375*t.^2+0.125*t+3;
plot(t,r,'-','LineWidth',2,'Color',[0 0 0])
t = linspace(0.8262,1.6639,100)';
r = 0.1333*t.^2-1.0333;
plot(t,r,'-','LineWidth',2,'Color',[0 0 0])
t = linspace(-0.9785,1.6639,100)';
r = 1-t;
plot(t,r,'-','LineWidth',2,'Color',[0 0 0])
t = linspace(-1.7476,0.8262,100)';
r = -1/3*t-2/3;
plot(t,r,'-','LineWidth',2,'Color',[0 0 0])
x(2)+0.9375*x(1).^2-0.125*x(1)-3;

---

A continuación, resolvemos el problema de optimización con la función "fmincon", como muestro a continuación:

A = [1 1; -1/3 -1];
B = [1; 2/3];
[xs,fs] = fmincon(@(u) f01(u), [0,0],A,B,[],[],[],[],@(u) frest01(u))

Matlab devuelve la solución:

xs =

    0.6512    0.3488


fs =

    0.1330

Finalmente, representamos la solución del problema en el gráfico:

plot(xs(1),xs(2),'o','MarkerSize',10,'Color',[1 0 0],'LineWidth',2)

---

Por último, vamos a añadir la siguiente restricción de igualdad no lineal.

• y - (0.08·x³ + 0.4·x² - 0.4·x - 1.2)·(x+2 )^-1 = 0

Es decir, la solución debe encontrarse, sobre dicha curva.

t = (-2:0.05:2)';
r = (0.08*t.^3+0.4*t.^2-0.4*t-1.2)./(t+2);
plot(t,r,'-','LineWidth',2,'Color',[1 0 0])

Para ello, lo único que tenemos que hacer es modificar el archivo con las restricciones no lineales ("frest01"), como muestro a continuación:

function [z zz] = frest01(x)

z(1) = x(2)+0.9375*x(1).^2-0.125*x(1)-3;

z(2) = -x(2)+0.1333*x(1).^2-1.0333;

zz = x(2) - (0.08*x(1).^3+0.4*x(1).^2-0.4*x(1)-1.2)./(x(1)+2);

Finalmente, calculamos la nueva solución:

[xs,fs] = fmincon(@(u) f01(u), [0,0],A,B,[],[],[],[],@(u) frest01(u))

plot(xs(1),xs(2),'o','MarkerSize',10,'Color',[0 0 0],'LineWidth',2)

Matlab calcula los siguientes valores:

xs =

    0.3428   -0.5493

fs =

    1.3212