%% Optimización con restricciones de igualdad y desigualdad.
% Copyright (c) 2022 G.T.A.S. (Grupo Tratamiento Avanzado de Señal). Universidad de Cantabria.
%
%% Índice:
%%     - Introducción y planteamiento del problema
%%     - Optimización con una restricción de igualdad y una de desigualdad
%                 - Definiciones
%                 - Condiciones de Karush-Kuhn-Tucker (KKT)
%                 - Clasificación de los puntos y puntos solución
%                 - Ejemplos
%%    - Optimización con una restricción de igualdad y varias restricciones de desigualdad
%                 - Definiciones
%                 - Condiciones de Karush-Kuhn-Tucker (KKT)
%                 - Clasificación de los puntos y puntos solución
%                 - Ejemplos
%%   - Métodos numéricos: función fmincon
%                 - Ejemplos
%%   - Ejercicios propuestos
%
%
%% Introducción y planteamiento del problema
%
% Determinar el mínimo de una función f(\mathbf x)   (función objetivo) donde la variable independiente \mathbf x está sujeta a restricciones de igualdad y de desigualdad: g_i(\mathbf x) =  0, \ \  i=1,\ldots,K_g  y   h_j(\mathbf x) \leq  0, \ \  j=1,\ldots,K_h
% Las funciones f(\mathbf x), \ \{g_i(\mathbf x)\}, \ \{h_j(\mathbf x) \} están definidas en un subconjunto abierto D en \mathbb{R}^M:    f(\mathbf x) :  D \subset \mathbb{R}^M \longrightarrow \mathbb{R}, \quad \quad
% g_i(\mathbf x) :  D \subset \mathbb{R}^M \longrightarrow \mathbb{R}, \quad \quad
% h_j(\mathbf x) :  D \subset \mathbb{R}^M \longrightarrow \mathbb{R},
% Se supone que todas las funciones f(\mathbf x), \ \{g_i(\mathbf x)\}, \ \{h_j(\mathbf x) \} son C^1, es decir, tienen derivadas parciales continuas en D.
%
% Notación en formato estandar:
%                 \min \quad f(\mathbf{x}) \\
% \matrix{s.t. \quad g_1(\mathbf x) = 0 \cr
%        \vdots \cr
% \quad \quad g_{Kg} (\mathbf x) = 0 \cr
% \quad \quad h_1(\mathbf x) \leq 0 \cr
% \vdots \cr
% \quad \quad h_{Kh}(\mathbf x) \leq 0
% }
%
%%
%% Optimización con una restricción de igualdad y una de desigualdad
%
%                 \min \quad f(\mathbf{x}) \\
% s.t. \quad g(\mathbf x) = 0 \\
% \quad \quad h(\mathbf x) \leq 0
%
%%  Definiciones:
%
% - Puntos factibles:   Son los puntos \mathbf x que cumplen las restricciones,  g(\mathbf x) = 0, \quad h(\mathbf{x}) \leq 0
%                 Se pueden clasificar en dos tipos teniendo en cuenta la restricción de desigualdad:
%                1)   Puntos que no saturan la restricción (ó puntos interiores):     h(\mathbf{x}) < 0
%                 2)   Puntos que saturan la restricción (ó puntos frontera):   h(\mathbf{x}) = 0
%
% - Puntos irregulares:  Son puntos factibles en los que los gradientes de las restricciones asociadas a las condiciones saturadas son linealmente dependientes. Dos posibilidades:
% 1) Puntos irregulares interiores:  h(\mathbf x) < 0, \ g(\mathbf x) = 0,  \ \nabla g(\mathbf x) = \mathbf 0
%  2) Puntos irregulares frontera :  h(\mathbf x) = 0,  g(\mathbf x) = 0,  \{\nabla g(\mathbf x),\nabla h(\mathbf x)\}
%   son linealmente dependientes.
%
% - Puntos KKT: Son puntos factibles que cumplen las llamadas condiciones de Karush-Kuhn-Tucker (KKT)
%
%% Condiciones necesarias de Karush-Kuhn-Tucker (KKT):
%
%                 \left\{
% \matrix{
% \nabla f(\mathbf{x}) +  \lambda \nabla g(\mathbf{x}) +  \mu \nabla h(\mathbf{x}) = \mathbf 0 \cr
% g(\mathbf x) = 0 \cr
% \mu \ h(\mathbf{x}) = 0 \cr
% \mu \geq 0, \cr
% h(\mathbf{x}) \leq 0
% }
%         donde \lambda, \mu son números reales que se conocen como multiplicador de Lagrange y multiplicador KKT, respectivamente.
% Además, de acuerdo a la clasificación de puntos factibles (puntos interiores y puntos frontera), las condiciones KKT se pueden particularizar para los puntos frontera y para los puntos interiores:
% 1) Puntos KKT interiores:  h(\mathbf x) < 0  \Rightarrow \quad  \mu = 0
%                          \Rightarrow
% \left\{
% \matrix{
% \nabla f(\mathbf{x})  + \lambda  \nabla g(\mathbf x) = \mathbf 0 \cr
% g(\mathbf{x}) = 0 \cr
% h(\mathbf{x}) < 0
% }
% 2) Puntos KKT frontera:\quad h(\mathbf x) = 0
%                          \Rightarrow
% \left\{
% \matrix{
% \nabla f(\mathbf{x}) +  \lambda \nabla g(\mathbf{x}) + \mu \nabla h(\mathbf x)= \mathbf 0 \cr
% g(\mathbf{x})= 0 \cr
% h(\mathbf{x})= 0 \cr
% \mu \geq 0
%
% }
%
%% Clasificación de los puntos:
%
% D:
% \left\{
% \matrix{
% - \ \textnormal{Puntos factibles}:
% \left\{\matrix{
% - \ \textnormal{Puntos regulares}:
% \left\{\matrix{
% - \ \textnormal{Puntos KKT} ^{(*)} \
% \cr
% \cr
% - \ \textnormal{Puntos no KKT}
% }
% \cr
% \cr
% - \ \textnormal{Puntos irregulares}^{(*)} \quad \quad \quad \quad \quad \quad \quad \quad
% }
% \cr
% \cr
% - \ \textnormal{Puntos no factibles} \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad
% }
%
% (*) Si el problema tiene solución, esta pertenece al conjunto de los puntos irregulares o al conjunto de los puntos KKT (interiores o frontera). Normalmente, estos conjuntos contienen un número finito de puntos, o ninguno.
% En algunos casos muy particulares puede haber puntos irregulares que también cumplan las condiciones KKT. De acuerdo a la clasificación anterior, nosotros los consideraremos como puntos irregulares.
% Las condiciones KKT son necesarias, pero no suficientes:  Si \mathbf x^* es solución del problema y es un punto regular, entonces es un punto KKT. Pero puede haber puntos KKT que no sean solución.
%
%
%% Ejemplo 1
%
%                     \min \quad x_1^2+x_2^2 \\
%  s.t. \quad x_1 + x_2 \geq 1\\
%  \quad \quad x_2 = 2 x_1
% Solución:
% D = \mathbb{R}^2.
% En primer lugar se escribe el problema en la forma estandar:
%                     \min \quad x_1^2+x_2^2 \\
%  s.t. \quad 1-x_1-x_2 \leq 0 \\
%   \quad \quad 2 x_1-x_2 = 0
% Por tanto,   f(\mathbf x) = x_1^2 + x_2^2, \quad \quad
% h(\mathbf x) = 1 - x_1 - x_2, \quad \quad
% g(\mathbf x) = 2 x_1-x_2
%
% Puntos irregulares:
% 1) Puntos irregulares interiores:
% \nabla g(\mathbf x) =\left[\matrix{
% 2 \cr
% -1
% }\right]. Por tanto, no hay puntos irregulares interiores.
% 2) Puntos irregulares frontera:
% \nabla h(\mathbf x) =\left[\matrix{
% -1 \cr
% -1
% }\right]  y  \nabla g(\mathbf x) =\left[\matrix{
% 2 \cr
% -1
% }\right] son linealmente independientes. Por tanto, no hay puntos irregulares frontera.
%
% Puntos KKT:
% 1) Puntos KKT interiores:
% \Rightarrow
% \left\{
% \matrix{
% \nabla f(\mathbf{x})  + \lambda \nabla g(\mathbf x) = \mathbf 0 \cr
% g(\mathbf{x}) = 0 \cr
% h(\mathbf{x}) < 0
% }     \Rightarrow \left\{ \matrix{
% 2 x_1 + 2 \lambda = 0 \cr
% 2 x_2 - \lambda  = 0 \cr
% 2 x_1 - x_2  = 0 \cr
% 1-x_1-x_2 < 0
% }
% Las tres primeras ecuaciones constituyen un sistema lineal,  \mathbf A \mathbf y = \mathbf b, donde
% \mathbf A = \left( \matrix{
% 2 & 0 & 2 \cr
% 0 & 2 & -1 \cr
% 2 & -1 & 0
% } \right ) ,  \quad \quad
% \mathbf y = \left( \matrix{x_1 \cr x_2 \cr \lambda} \right ) ,  \quad \quad
% \mathbf b = \left( \matrix{0 \cr 0 \cr 0} \right )

% Resuelve el sistema de ecuaciones
A=[2 0 2;0 2 -1; 2 -1 0]; b=[0 0 0]';
disp(['|A| = ',num2str(det(A))])
y = A\b;
disp(['solución: ',num2str(y')])

% Pero este punto no es un punto interior a la restricción de desigualdad (no cumple la última condición). Por tanto, no hay puntos KKT interiores.
%
% 2) Puntos KKT frontera:
%  \left\{
% \matrix{
% \nabla f(\mathbf{x}) +  \lambda \nabla g(\mathbf{x}) + \mu \nabla h(\mathbf x)= \mathbf 0 \cr
% g(\mathbf{x})= 0 \cr
% h(\mathbf{x})= 0 \cr
% \mu \geq 0
% }  \Rightarrow      \left\{ \matrix{
% 2 x_1 + 2 \lambda - \mu = 0 \cr
% 2 x_2 - \lambda - \mu = 0 \cr
% 2 x_1 - x_2  = 0 \cr
% 1-x_1-x_2 = 0 \cr
% \mu \geq 0
% }
% Las cuatro primeras ecuaciones constituyen un sistema lineal,  \mathbf A \mathbf y = \mathbf b, donde
% \mathbf A = \left( \matrix{
% 2 & 0 & 2  &-1 \cr
% 0 & 2 & -1 & -1\cr
% 2 & -1 & 0 & 0 \cr
% -1 & -1 & 0 & 0
% } \right ) ,  \quad \quad
% \mathbf y = \left( \matrix{x_1 \cr x_2 \cr \lambda \cr \mu} \right ) ,  \quad \quad
% \mathbf b = \left( \matrix{0 \cr 0 \cr 0 \cr -1} \right )
% Como \det(\mathbf A) \neq 0, el sistema tiene una única solución:

% Resuelve el sistema de ecuaciones
A=[2 0 2 -1; 0 2 -1 -1; 2 -1 0 0; -1 -1 0 0]; b=[0 0 0 -1]';
y = A\b;
disp(['solución: ',num2str(y')])

% Este punto cumple la última condición pòr lo que es un punto KKT frontera.
% Gráficamente se ve que es la solución del problema.

funcion = @(x1,x2) x1.^2 + x2.^2;
restriccion_h = @(x1,x2) 1 - x1 - x2;
restriccion_g = @(x1,x2) 2*x1 - x2;
region = [-1 1 -1 1];

fcontour(funcion, region, 'LevelStep', 0.04)
xlabel('x_1'); ylabel('x_2')
colorbar
hold on
fimplicit(restriccion_h, region,'r')
fimplicit(restriccion_g, region,'k')
plot(y(1),y(2),'ko')
title('f(x), g(x)=0, h(x)=0')
hold off

%
%
%% Ejemplo 2
%
% Considere la familia de triangulos rectángulos de area A=2 en los que un cateto es mayor del doble que el otro. Determine cuál de ellos tiene menor hipotenusa.
%
% Solución:
% D = \mathbb{R}^2.
%                                 \min \quad x_1^2 +  x_2^2 \\
%  s.t. \quad x_1  x_2 -2 A=0 \\
% \quad \quad -x_1 + 2 x_2 \leq 0
%
% En este caso   f(\mathbf x) =  x_1^2 + x_2^2,
% \quad  g(\mathbf x)= x_1 x_2 - 2A,
% \quad h(\mathbf x)=-x_1 + 2 x_2
%
% Puntos irregulares:
% 1) Puntos irregulares interiores:
% \nabla g(\mathbf x) =\left[\matrix{
% x_2 \cr
% x_1
% }\right] = \mathbf 0 en \mathbf x =\mathbf 0, pero no es un punto factible.
% 1) Puntos irregulares frontera:
% h(\mathbf x) =0  y  g(\mathbf x) =0 solo se cumple en \mathbf x =\sqrt A
% \left[\matrix{
% 2 \cr
% 1
% }\right]. Los gradientes en ese punto son
% \nabla h(\mathbf x) =\left[\matrix{
% -1 \cr
% 2
% }\right],  \nabla g(\mathbf x) = \sqrt A \left[\matrix{
% 1 \cr
% 2
% }\right], que son linealmente independientes.

rank([-1,1;2 2])

% Por tanto, no hay puntos irregulares.
%
% Puntos KKT:
% 1) Puntos KKT interiores:
%   \Rightarrow
% \left\{
% \matrix{
% \nabla f(\mathbf{x})  + \lambda  \nabla g(\mathbf x) = \mathbf 0 \cr
% g(\mathbf{x}) = 0 \cr
% h(\mathbf{x}) < 0
% }     \Rightarrow \left\{ \matrix{
% 2 x_1 + \lambda  x_2 = 0 \cr
% 2 x_2 + \lambda  x_1 = 0 \cr
% x_1 x_2  -2A = 0 \cr
% -x_1 + 2 x_2< 0
% }
% Las 3 primeras condiciones constituyen un sistema de ecuaciones no lineales. Resuelvo

A=2;
sistema = @(y)[2*y(1)+y(2)*y(3) , 2*y(2)+y(1)*y(3) , y(1)*y(2)-2*A];
y0=[1,1,0];
y = fsolve(sistema,y0);
disp(['solución: ',num2str(y)])
disp(['valor del sistema de ecuaciones en la solución: ',num2str(sistema(y))])

%  Pero la solución, \mathbf x = (2,2), no cumple la última condición.
%
% 2) Puntos KKT frontera:
%
% \left\{
% \matrix{
% \nabla f(\mathbf{x}) +  \lambda \nabla g(\mathbf{x}) + \mu \nabla h(\mathbf x)= \mathbf 0 \cr
% g(\mathbf{x})= 0 \cr
% h(\mathbf{x})= 0 \cr
% \mu \geq 0
% }  \Rightarrow      \left\{ \matrix{
% 2 x_1 + \lambda x_2 - \mu = 0 \cr
% 2 x_2 + \lambda x_1 + 2 \mu = 0 \cr
% x_1 x_2 -2A = 0 \cr
% -x_1 + 2x_2 = 0 \cr
% \mu \geq 0
% }
% Las 3 primeras condiciones constituyen un sistema se ecuaciones no lineales:

sistema = @(y)[2*y(1)+y(2)*y(3)-y(4), 2*y(2)+y(1)*y(3)+2*y(4), y(1)*y(2)-2*A, -y(1)+2*y(2)];
y0=[1,1,0,0];
y = fsolve(sistema,y0);
disp(['solución: ',num2str(y)])
disp(['valor del sistema de ecuaciones en la solución: ',num2str(sistema(y))])

%  La solución, \mathbf x = (2.8284,1.4142) con \mu = 2.1213 \geq 0  cumple la última condición, por tanto es un punto KKT frontera. Se comprueba gráficamente que es la solución:

funcion = @(x1,x2) x1.^2 + x2.^2;
restriccion1 = @(x1,x2) x1.*x2-2*A ;
restriccion2 = @(x1,x2) -x1+2*x2 ;
region = [0 4 0 4];

fcontour(funcion, region, 'LevelStep', 0.2)
colorbar
xlabel('x_1'); ylabel('x_2')

hold on
fimplicit(restriccion1, region,'k')
fimplicit(restriccion2, region,'r')
title('f(x), g(x)=0, h(x)=0')
plot(y(1),y(2),'ko')
hold off

%
%
%% Ejemplo 3
%
%                     \min \quad 2 x_1 + 2 x_2 + 2 x_3 \\
%  s.t. \quad x_1^2 + x_2^2 -x_3 -1 \leq 0\\
%  \quad \quad (x_1-3)^2 + x_2^2 + 2 x_3 - 4 = 0
% Solución:
% D = \mathbb{R}^3.
% f(\mathbf x) = x_1+x_2+x_3, \quad \quad
% h(\mathbf x) = x_1^2 + x_2^2 -x_3 -1, \quad \quad
% g(\mathbf x) = (x_1-3)^2 + x_2^2 + 2 x_3 - 4
%
% Puntos irregulares:
% 1) Puntos irregulares interiores:
% \nabla g(\mathbf x) =\left[\matrix{
% 2 (x_1-3) \cr
% 2 x_2 \cr
% 2
% }\right]
% \neq \mathbf 0. Por tanto no hay puntos irregulares interiores.
%
% 1) Puntos irregulares frontera:
% Los gradientes deben ser linealmente dependientes \nabla h(\mathbf x) = \alpha \nabla g(\mathbf x)  \Rightarrow
% \left[\matrix{
% 2 x_1 \cr
% 2 x_2 \cr
% -1
% }\right] = \alpha
% \left[\matrix{
% 2 (x_1-3) \cr
% 2 x_2 \cr
% 2
% }\right]
% Este es un sistema (no lineal) de 3 ecuaciones con 3 incognitas, x_1,x_2,\alpha, que se puede resolver analíticamente: De la última ecuación se tiene \alpha = -1/2. De la segunda ecuación se tiene x_2 = 0. Y de la primera se tiene x_1 = 1.
% Además, se debe cumplir  h(\mathbf x) =0  y  g(\mathbf x) =0. Para que se cumpla h(\mathbf x) =0, se tiene que cumplir que x_3=0, y también se cumple g(\mathbf x) =0. Por tanto \mathbf x = (1,0,0)es un punto irregular frontera. En este punto la función de coste es f(1,0,0) = 1.
%
% Puntos KKT:
% 1) Puntos KKT interiores:
%   \Rightarrow
% \left\{
% \matrix{
% \nabla f(\mathbf{x})  + \lambda  \nabla g(\mathbf x) = \mathbf 0 \cr
% g(\mathbf{x}) = 0 \cr
% h(\mathbf{x}) < 0
% }     \Rightarrow \left\{ \matrix{
% 2 + \lambda 2 (x_1-3) = 0 \cr
% 2 + \lambda  2 x_2 = 0 \cr
% 2 + 2 \lambda  = 0 \cr
% (x_1-3)^2 + x_2^2 + 2 x_3-4  = 0 \cr
% x_1^2 + x_2^2 - x_3 -1< 0
% }
% Las 3 primeras condiciones constituyen un sistema de ecuaciones no lineales que se puede resolver analiticamente:  \lambda = -1, \ x_2=1, \ x_1 = 4. Estonces, por la cuarta ecuación, x_3 = 1. Pero entonces no se cumple la condición de desigualdad. Por tanto, no hay puntos KKT interiores.
%
% 2) Puntos KKT frontera:
%
% \left\{
% \matrix{
% \nabla f(\mathbf{x}) +  \lambda \nabla g(\mathbf{x}) + \mu \nabla h(\mathbf x)= \mathbf 0 \cr
% g(\mathbf{x})= 0 \cr
% h(\mathbf{x})= 0 \cr
% \mu \geq 0
% }  \Rightarrow      \left\{ \matrix{
% 2+  \lambda 2 (x_1-3) + \mu 2 x_1 = 0 \cr
% 2 + \lambda 2 x_2 + \mu 2 x_2 = 0 \cr
% 2 + \lambda 2  - \mu = 0 \cr
% (x_1-3)^2 + x_2^2 + 2 x_3-4  = 0 \cr
% x_1^2 +  x_2^2  -x_3 -1 = 0 \cr
% \mu \geq 0
% }
%

sistema = @(y)[2 + y(4)*2*(y(1)-3) + y(5)*2*y(1),...
               2 + y(4)*2*y(2) + y(5)*2*y(2),...
               2 + y(4)*2-y(5),...
               y(1)^2+y(2)^2-y(3)-1,...
               (y(1)-3)^2+y(2)^2+2*y(3)-4];
y0=[0,0,0,0,0];
y = fsolve(sistema,y0);
disp(['solución: ',num2str(y)])
disp(['valor del sistema de ecuaciones en la solución: ',num2str(sistema(y))])

% No converge a una solución despues de probar con distintos puntos semilla. El algoritmo se detiene cerca al punto irregular \mathbf x = (1,0,0), pero en este punto no se cumple la segunda ecuación. Por tanto no hay puntos KKT.
% La solución es el punto irregular \mathbf x = (1,0,0).
%
%
%
%% Optimización con una restricción de igualdad y varias restricciones de desigualdad
%
%                 \min \quad f(\mathbf{x}) \\
% s.t. \quad g(\mathbf x) = 0 \\
% \quad \quad h_1(\mathbf x) \leq 0 \\
% \quad \quad \vdots \\
% \quad \quad h_{K}(\mathbf x) \leq 0
%
%% Definiciones:
%
% - Puntos factibles:   Son los puntos que cumplen las restricciones,  g(\mathbf x) = 0, \quad h_1(\mathbf{x}) \leq 0, \  \ldots, \quad h_K(\mathbf{x}) \leq 0
%                 Para cada restricción de desigualdad, los puntos factibles se pueden clasificar en dos tipos:
%                 -   Puntos que saturan la restricción (ó puntos frontera):      h_k(\mathbf{x}) = 0
%                 -   Puntos que no saturan la restricción (ó puntos interiores):     h_k(\mathbf{x}) < 0
%
%                 Entonces los puntos factibles se pueden clasificar en distintos grupos atendiendo a cuales son las restricciones que saturan.
%                 Por ejemplo, si el problema tiene K=2 restricciones de desigualdad, los puntos interiores se pueden clasificar en 2^K = 4 tipos:
%                 1) Puntos factibles interiores para ambas restricciones:  g(\mathbf x) = 0, \ h_1(\mathbf{x}) < 0, \ h_2(\mathbf{x}) < 0
%                 2) Puntos factibles frontera para ambas restricciones:  g(\mathbf x) = 0, \ h_1(\mathbf{x}) = 0, \ h_2(\mathbf{x}) = 0
%                 3) Puntos factibles frontera para la primera restricción e interior para la segunda:  g(\mathbf x) = 0, \ h_1(\mathbf{x}) = 0, \ h_2(\mathbf{x}) < 0
%                 4) Puntos factibles interiores para la primera restricción y frontera para la segunda:  g(\mathbf x) = 0, \ h_1(\mathbf{x}) < 0, \ h_2(\mathbf{x}) = 0
%
% - Puntos irregulares:  Son puntos factibles en los que los gradientes de las restricciones que satura (incluyendo el gradiente de la restricción de igualdad) son vectores linealmente dependientes.
%                 Por ejemplo, si el problema tiene K=2 restricciones de desigualdad, los puntos irregulares pueden ser de 2^K = 4 tipos:
%                 1) Puntos irregulares interiores para ambas restricciones:  g(\mathbf{x}) = 0, \ h_1(\mathbf{x}) < 0, \ h_2(\mathbf{x}) < 0, \ \nabla g(\mathbf{x})  = \mathbf 0
%                 2) Puntos irregulares frontera para ambas restricciones:  g(\mathbf{x}) = 0, \ h_1(\mathbf{x}) = 0, \ h_2(\mathbf{x}) = 0, \
% \{ \nabla g(\mathbf{x}),   \nabla h_1(\mathbf{x}), \nabla h_2(\mathbf{x}) \}  son vectores linealmente dependientes.
%                 3) Puntos irregulares frontera para la primera restricción e interior para la segunda:  g(\mathbf{x}) = 0, \ h_1(\mathbf{x}) = 0, \ h_2(\mathbf{x}) < 0, \
% \{ \nabla g(\mathbf{x}),   \nabla h_1(\mathbf{x}) \}  son vectores linealmente dependientes.
%                 4) Puntos irregulares interiores para la primera restricción y frontera para la segunda:  g(\mathbf{x}) = 0, \ h_1(\mathbf{x}) < 0, \ h_2(\mathbf{x}) = 0, \
% \{ \nabla g(\mathbf{x}),   \nabla h_2(\mathbf{x}) \}  son vectores linealmente dependientes
%
% - Puntos KKT: Son puntos factibles que cumplen las llamadas condiciones de Karush-Kuhn-Tucker (KKT).
%
%% Condiciones necesarias de Karush-Kuhn-Tucker (KKT):
%
%                 \left\{
% \matrix{
% \nabla f(\mathbf{x}) +  \lambda \nabla g(\mathbf{x}) +  \sum_{k=1}^K \mu_k \nabla h_k(\mathbf{x}) = \mathbf 0 \cr
% g(\mathbf x) = 0 \cr
% \mu_1 h_1(\mathbf{x}) = 0 \cr
% \mu_1 \geq 0, \cr
% h_1(\mathbf{x}) \leq 0 \cr
% \vdots \cr
% \mu_K h_K(\mathbf{x}) = 0 \cr
% \mu_K \geq 0, \cr
% h_K(\mathbf{x}) \leq 0
% }
%
%         donde (\lambda, \mu_1, \ldots, \mu_K) son números reales. El primero se conoce como multiplicador de Lagrange y el resto como multiplicadores KKT.
% Las condiciones KKT se pueden particularizar en función de que los puntos sean interiores o frontera para cada restricción de desigualdad. Por ejemplo, en el caso de K=2 restricciones de desigualdad:
%
% 1) Puntos KKT interiores a ambas restricciones:  h_1(\mathbf x) <0,  \ h_2(\mathbf x) <0  \Rightarrow \quad  \mu_1 = 0, \ \mu_2=0
%                          \Rightarrow
% \left\{
% \matrix{
% \nabla f(\mathbf{x})  + \lambda  \nabla g(\mathbf x) = \mathbf 0 \cr
% g(\mathbf{x}) = 0 \cr
% h_1(\mathbf{x}) < 0 \cr
% h_2(\mathbf{x}) < 0
% }
% 2) Puntos KKT frontera para ambas restricciones:    h_1(\mathbf x)=0, \ \ h_2(\mathbf x) = 0
%                          \Rightarrow
% \left\{
% \matrix{
% \nabla f(\mathbf{x})  + \lambda \nabla g(\mathbf x) + \mu_1 \nabla h_1(\mathbf x) + \mu_2 \nabla h_2(\mathbf x) = \mathbf 0 \cr
% g(\mathbf{x}) = 0 \cr
% h_1(\mathbf{x}) = 0 \cr
% h_2(\mathbf{x}) = 0 \cr
% \mu_1 \geq 0 \cr
% \mu_2 \geq 0
% }
% 3) Puntos KKT frontera para la primera restricción e interiores para la segunda:   h_1(\mathbf x) = 0, \ \ h_2(\mathbf x) < 0 \quad  \Rightarrow \mu_2=0
%                                    \Rightarrow
% \left\{
% \matrix{
% \nabla f(\mathbf{x})  + \lambda \nabla g(\mathbf x)+ \mu_1 \nabla h_1(\mathbf x)= \mathbf 0 \cr
% g(\mathbf{x}) = 0 \cr
% h_1(\mathbf{x}) = 0 \cr
% h_2(\mathbf{x}) < 0 \cr
% \mu_1 \geq 0
% }
% 4) Puntos KKT interiores para la primera restricción y frontera para la segunda:   h_1(\mathbf x) < 0, \ \ h_2(\mathbf x) = 0 \quad  \Rightarrow \mu_1=0
%                                    \Rightarrow
% \left\{
% \matrix{
% \nabla f(\mathbf{x})  + + \lambda \nabla g(\mathbf x) + \mu_2 \nabla h_2(\mathbf x)= \mathbf 0 \cr
% g(\mathbf{x}) = 0 \cr
% h_1(\mathbf{x}) < 0 \cr
% h_2(\mathbf{x}) = 0 \cr
% \mu_2 \geq 0
% }
%
% En general, si hubiera Krestricciones de desigualdad, habría 2^K tipos posibles de puntos KKT.
%
%% Clasificación de los puntos:
%
% D:
% \left\{
% \matrix{
% - \ \textnormal{Puntos factibles}:
% \left\{\matrix{
% - \ \textnormal{Puntos regulares}:
% \left\{\matrix{
% - \ \textnormal{Puntos KKT} ^{(*)} \
% \cr
% - \ \textnormal{Puntos no KKT}
% }
% \cr
% \cr
% - \ \textnormal{Puntos irregulares}^{(*)} \quad \quad \quad \quad \quad \quad \quad
% }
% \cr
% \cr
% - \ \textnormal{Puntos no factibles} \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad
% }
%
% (*) Si el problema tiene solución, esta pertenece al conjunto de los puntos irregulares o al conjunto de los puntos KKT (interiores o frontera). Normalmente, estos conjuntos contienen un número finito de puntos, o ninguno.
% En algunos casos muy particulares puede haber puntos irregulares que también cumplan las condiciones KKT. De acuerdo a la clasificación anterior, estos se considerarán irregulares.
% Las condiciones KKT son necesarias, pero no suficientes:  Si \mathbf x^* es solución del problema y es un punto regular, entonces es un punto KKT. Pero puede haber puntos KKT que no sean solución.
%%
%
%% Ejemplo 4
%
%                     \min \quad -x_1^2 + 2 x_1 + 4 x_2 \\
%  s.t. \quad -x_1 +4 x_2  = 0 \\
% \quad \quad x_1+x_2-4 \leq 0 \\
% \quad \quad -x_1+2 \leq 0
% Solución:
% D = \mathbb{R}^2.
% Ya está en forma estandar.
% Los gradiantes de las funciones involucradas:
% \nabla g(\mathbf{x})  = \left( \matrix{-1 \cr 4} \right), \quad
% \nabla h_1(\mathbf{x})  = \left( \matrix{1 \cr 1} \right), \quad
% \nabla h_2(\mathbf{x})  = \left( \matrix{-1 \cr 0} \right), \quad
% \nabla f(\mathbf{x}) = \left( \matrix{- 2 x_1 +2 \cr 4} \right)
%
% Puntos irregulares:
% 1)   g(\mathbf x) = 0, \ h_1(\mathbf x) < 0, h_2(\mathbf x) < 0, \nabla g(\mathbf x) = \mathbf 0
% No se cumple la última condición. Por tanto, no hay puntos irregulares interiores a ambas restricciones de desigualdad.
%
% 2)   g(\mathbf{x}) = 0, \ h_1(\mathbf{x}) = 0, \ h_2(\mathbf{x}) = 0, \
% \{ \nabla g(\mathbf{x}),   \nabla h_1(\mathbf{x}), \nabla h_2(\mathbf{x}) \}  lin.  dep..
% No hay puntos factibles de este tipo. Por tanto, no hay puntos irregulares que saturan ambas restricciones de desigualdad.
%
% 3)    g(\mathbf{x}) = 0, \ h_1(\mathbf{x}) = 0, \ h_2(\mathbf{x}) < 0, \
% \{ \nabla g(\mathbf{x}),   \nabla h_1(\mathbf{x}) \}  lin. dep.
% \nabla h_1(\mathbf x) , \
% \nabla g(\mathbf x)   son linealmente independientes. Por tanto no hay puntos irregulares entre los puntos que son frontera para la primera restricción e interiores para la segunda.
%
% 4)   g(\mathbf{x}) = 0, \ h_1(\mathbf{x}) < 0, \ h_2(\mathbf{x}) = 0, \
% \{ \nabla g(\mathbf{x}),   \nabla h_2(\mathbf{x}) \}  lin. dep.
% \nabla h_2(\mathbf x), \nabla g(\mathbf x)  son linealmente independientes. Por tanto no hay puntos irregulares entre los puntos que son interiores para la primera restricción y frontera para la segunda.
%
%
%
% Puntos KKT:
% 1) h_1(\mathbf x) <0,  \ h_2(\mathbf x) <0  \Rightarrow \quad  \mu_1 = 0, \ \mu_2=0
%                          \left\{
% \matrix{
% \nabla f(\mathbf{x})  + \lambda  \nabla g(\mathbf x) = \mathbf 0 \cr
% g(\mathbf{x}) = 0 \cr
% h_1(\mathbf{x}) < 0 \cr
% h_2(\mathbf{x}) < 0
% } \Rightarrow
% \left\{
% \matrix{
% -2 x_1 +2- \lambda = 0 \cr
% 4 + 4 \lambda = 0 \cr
% -x_1+ 4 x_2 = 0 \cr
% x_1+x_2-4 \leq 0 \cr
% -x_1+2 \leq 0
% }
% De las tres primeras ecuaciones, \lambda = -1, \  x_1 = 3/2, \ x_2 = 3/8. Pero no se cumple la última desigualdad. No hay puntos KKT interiores a ambas restricciones.
% 2)   h_1(\mathbf x)=0, \ \ h_2(\mathbf x) = 0
%                         \Rightarrow
% \left\{
% \matrix{
% \nabla f(\mathbf{x})  + \lambda \nabla g(\mathbf x) + \mu_1 \nabla h_1(\mathbf x) + \mu_2 \nabla h_2(\mathbf x) = \mathbf 0 \cr
% g(\mathbf{x}) = 0 \cr
% h_1(\mathbf{x}) = 0 \cr
% h_2(\mathbf{x}) = 0 \cr
% \mu_1 \geq 0 \cr
% \mu_2 \geq 0
% }           \Rightarrow
% \left\{
% \matrix{
% -2 x_1 +2 -\lambda + \mu_1 - \mu_2 = 0 \cr
% 4+4\lambda + \mu_1 = 0 \cr
% -x_1+4x_2 = 0 \cr
% x_1+x_2-4  = 0 \cr
% -x_1+2 = 0 \cr
% \mu_1 \geq 0 \cr
% \mu_2 \geq 0
% }
% Las 3 últimas ecuaciones no se pueden cumplir simultaneamente por lo que no existen puntos frontera (para ambas restricciones) que cumplan las condiciones KKT.
% 3)   h_1(\mathbf x) = 0, \ \ h_2(\mathbf x) < 0 \quad  \Rightarrow \mu_2=0
%                                    \Rightarrow
% \left\{
% \matrix{
% \nabla f(\mathbf{x})  + \lambda \nabla g(\mathbf x)+ \mu_1 \nabla h_1(\mathbf x)= \mathbf 0 \cr
% g(\mathbf{x}) = 0 \cr
% h_1(\mathbf{x}) = 0 \cr
% h_2(\mathbf{x}) < 0 \cr
% \mu_1 \geq 0
% }             \Rightarrow
% \left\{
% \matrix{
% -2 x_1 +2 -\lambda + \mu_1 = 0 \cr
% 4+4\lambda + \mu_1 = 0 \cr
% -x_1+4x_2 = 0 \cr
% x_1+x_2-4  = 0 \cr
% -x_1+2 < 0 \cr
% \mu_1 \geq 0
% }
% Las 4 ecuaciones constituyen un sistema lineal:

A=[-2 0 -1 1; 0 0 4 1; -1 4 0 0; 1 1 0 0]; b=[-2 -4 0 4]';
y = A\b;
disp(['solución: ',num2str(y')])

% La solución es el punto \mathbf x = (16/5, 4/5) con \mu_1=68/25. Por tanto, este es un punto KKT y, por tanto, es un candidato a ser solución del problema.
% 4)    h_1(\mathbf x) < 0, \ \ h_2(\mathbf x) = 0 \quad  \Rightarrow \mu_1=0
%                                    \Rightarrow
% \left\{
% \matrix{
% \nabla f(\mathbf{x})  + + \lambda \nabla g(\mathbf x) + \mu_2 \nabla h_2(\mathbf x)= \mathbf 0 \cr
% g(\mathbf{x}) = 0 \cr
% h_1(\mathbf{x}) < 0 \cr
% h_2(\mathbf{x}) = 0 \cr
% \mu_2 \geq 0
% }     \Rightarrow
% \left\{
% \matrix{
% -2 x_1 +2 -\lambda  - \mu_2 = 0 \cr
% 4+4\lambda = 0 \cr
% -x_1+4x_2 = 0 \cr
% -x_1+2 = 0 \cr
% x_1+x_2-4  < 0 \cr
% \mu_2 \geq 0
% }
% Las 4 ecuaciones constituyen un sistema lineal:

A=[-2 0 -1 -1; 0 0 4 0; -1 4 0 0; -1 0 0 0]; b=[-2 -4 0 -2]';
y = A\b;
disp(['solución: ',num2str(y')])

% La solución  del sistema de ecuaciones es el punto \mathbf x = (2, 1/2) con \mu_2=-1. Por tanto, no es un punto KKT.
%
% El único punto que cumple las condiciones KKT es \mathbf x = (16/5, 4/5), que es frontera para la primera restricción de desigualdad e interior para la segunda.
% Gráficamente se ve que es solución del problema:

funcion = @(x1,x2) -x1.^2+4*x2+2*x1;
restriccion1 = @(x1,x2) -x1+4*x2;
restriccion2 = @(x1,x2) x2+x1-4;
restriccion3 = @(x1,x2) -x1+2;
region = [0 4 0 4];

fcontour(funcion, region, 'LevelStep', 0.1)
colorbar
xlabel('x_1'); ylabel('x_2')
hold on
fimplicit(restriccion1, region,'k')
fimplicit(restriccion2, region,'r')
fimplicit(restriccion3, region,'r')
plot(16/5, 4/5, 'ko')
hold off
title('f(x), g(x)=0, h_1(x), h_2(x)')

%
%
%% Ejemplo 5
%
% La capacidad de un canal de comunicaciones es una medida de la velocidad máxima (en bits/s/Hz) a la que se puede transmitir a través de él, sin errores.
% El problema "waterfilling" es un problema clásico en sistemas de comunicaciones y teoría de la información. En estos sistemas se dispone de varios canales de comunicación y una potencia total disponible P. La capacidad del canal m (en bits/s/Hz) es \log(1+a_m x_m), donde a_m es la ganancia del canal m y x_m es la potencia con la que se transmite por el canal m. El objetivo es repartir la potencia total disponible, P, entre los canales para maximizar la suma de la velocidad máxima de transmisión por los canales:
% En el caso de dos canales  (m=2) el problema se puede formular así:
%                     \max \quad \log(1+a_1 x_1) + \log(1+a_2 x_2) \\
%  s.t. \quad x_1+x_2 = P \\
% \quad \quad x_1 \geq 0 \\
% \quad \quad x_2 \geq 0
% donde P, \ a_1 y a_2 son números reales positivos.
%
% Solución:
% D = \mathbb{R}^2.
% Escribo el problema en la forma estandar
%                     \min \quad -\log(1+a_1 x_1) - \log(1+a_2 x_2) \\
%  s.t. \quad x_1+x_2 -P= 0 \\
% \quad \quad -x_1 \leq 0 \\
% \quad \quad -x_2 \leq 0
% Los gradientes de las funciones:
% \nabla g(\mathbf{x})  = \left( \matrix{1 \cr 1} \right), \quad
% \nabla h_1(\mathbf{x})  = \left( \matrix{-1 \cr 0} \right), \quad
% \nabla h_2(\mathbf{x})  = \left( \matrix{0 \cr -1} \right), \quad
% \nabla f(\mathbf{x}) = \left( \matrix{\frac{-a_1}{1+a_1 x_1} \cr \frac{-a_2}{1+a_2 x_2}} \right)
% Puntos irregulares:
% 1)   h_1(\mathbf x) = 0, \ h_2(\mathbf x) < 0
% \nabla h_1(\mathbf x)
%   y \nabla g(\mathbf x)son linealmente independientes. Por tanto no hay puntos irregulares entre los puntos que son frontera para la primera restricción e interiores para la segunda.
% 2)   h_1(\mathbf x) < 0, \ h_2(\mathbf x) = 0
% \nabla h_2(\mathbf x) , \nabla g(\mathbf x) son linealmente independientes. Por tanto no hay puntos irregulares entre los puntos que son interiores para la primera restricción y frontera para la segunda.
% 3)   h_1(\mathbf x) = 0, h_2(\mathbf x) = 0
% Se cumple en \mathbf x = (0,0) ,  pero este punto no es factible ya que no cumple g(\mathbf x) = \mathbf 0.
%  Por tanto, no hay puntos irregulares que saturan ambas restricciones de desigualdad.
%
% Puntos KKT:
% 1)   h_1(\mathbf x) <0,  \ h_2(\mathbf x) <0  \Rightarrow \quad  \mu_1 = 0, \ \mu_2=0
%                          \Rightarrow
% \left\{
% \matrix{
% \nabla f(\mathbf{x})  + \lambda  \nabla g(\mathbf x) = \mathbf 0 \cr
% g(\mathbf{x}) = 0 \cr
% h_1(\mathbf{x}) < 0 \cr
% h_2(\mathbf{x}) < 0
% } \Rightarrow
% \left\{
% \matrix{
% \frac{-a_1}{1+a_1 x_1}+\lambda  = 0 \cr
% \frac{-a_2}{1+a_2 x_2}+\lambda  = 0 \cr
% x_1+x_2-P = 0 \cr
% -x_1  < 0 \cr
% -x_2 < 0
% }
%
% De las primeras 3 ecuaciones se obtiene analíticamente   x_1 =\frac{P}{2}-\frac{1}{2 a_1}+\frac{1} {2a_2}, \quad
% x_2 =\frac{P}{2}+\frac{1}{2 a_1}-\frac{1} {2a_2} .
% Y para que el punto anterior sea KKT debe cumplir las dos últimas condiciones:  P \geq \frac{1}{a_1}-\frac{1}{a_2}, \quad P \geq \frac{1}{a_2}-\frac{1}{a_1}.
%
% 2)    h_1(\mathbf x)=0, \ \ h_2(\mathbf x) = 0
%                         \Rightarrow
% \left\{
% \matrix{
% \nabla f(\mathbf{x})  + \lambda \nabla g(\mathbf x) + \mu_1 \nabla h_1(\mathbf x) + \mu_2 \nabla h_2(\mathbf x) = \mathbf 0 \cr
% g(\mathbf{x}) = 0 \cr
% h_1(\mathbf{x}) = 0 \cr
% h_2(\mathbf{x}) = 0 \cr
% \mu_1 \geq 0 \cr
% \mu_2 \geq 0
% }           \Rightarrow
% \left\{
% \matrix{
% \frac{-a_1}{1+a_1 x_1}+\lambda - \mu_1 = 0 \cr
% \frac{-a_2}{1+a_2 x_2}+\lambda  - \mu_2 = 0 \cr
% x_1+x_2-P = 0 \cr
% x_1  = 0 \cr
% x_2 = 0 \cr
% \mu_1 \geq 0 \cr
% \mu_2 \geq 0 \cr
% }
% De las condiciones 3-5 se ve que no hay ningún punto que cumpla las condiciones.
%
% 3)    h_1(\mathbf x) = 0, \ \ h_2(\mathbf x) < 0 \quad  \Rightarrow \mu_2=0
%                                    \Rightarrow
% \left\{
% \matrix{
% \nabla f(\mathbf{x})  + \lambda \nabla g(\mathbf x)+ \mu_1 \nabla h_1(\mathbf x)= \mathbf 0 \cr
% g(\mathbf{x}) = 0 \cr
% h_1(\mathbf{x}) = 0 \cr
% h_2(\mathbf{x}) < 0 \cr
% \mu_1 \geq 0
% } \Rightarrow
% \left\{
% \matrix{
% \frac{-a_1}{1+a_1 x_1}+\lambda - \mu_1 = 0 \cr
% \frac{-a_2}{1+a_2 x_2}+\lambda   = 0 \cr
% x_1+x_2-P = 0 \cr
% x_1  = 0 \cr
% -x_2 < 0 \cr
% \mu_1 \geq 0 \cr
% }
% Ver que de las ecuaciones 3-4 se tiene que \mathbf x =(0,P). Entonces, de las dos primeras ecuaciones obtenemos \mu_1 = -a_1+\frac{a_2}{1+a_2 P},  con lo que  debemos comprobar si esto cumple la última condición, es decir, si a_1\leq\frac{a_2}{1+a_2 P} \ \Rightarrow \ P \leq \frac{1}{a_1}-\frac{1}{a_2}.
% Por tanto, si P \leq \frac{1}{a_1}-\frac{1}{a_2}, el punto  \mathbf x =(0,P) es KKT.
%
% 4)    h_1(\mathbf x) < 0, \ \ h_2(\mathbf x) = 0 \quad  \Rightarrow \mu_1=0
%                                    \Rightarrow
% \left\{
% \matrix{
% \nabla f(\mathbf{x})  + + \lambda \nabla g(\mathbf x) + \mu_2 \nabla h_2(\mathbf x)= \mathbf 0 \cr
% g(\mathbf{x}) = 0 \cr
% h_1(\mathbf{x}) < 0 \cr
% h_2(\mathbf{x}) = 0 \cr
% \mu_2 \geq 0
% } \Rightarrow
% \left\{
% \matrix{
% \frac{-a_1}{1+a_1 x_1}+\lambda  = 0 \cr
% \frac{-a_2}{1+a_2 x_2}+\lambda - \mu_2  = 0 \cr
% x_1+x_2-P = 0 \cr
% -x_1  < 0 \cr
% x_2 = 0 \cr
% \mu_2 \geq 0 \cr
% }
% Ver que de las condiciones 3 y 5 se tiene que \mathbf x =(P,0). Entonces, de las dos primeras ecuaciones se obtiene \mu_2 = -a_2+\frac{a_1}{1+a_1 P},  con lo que  debemos comprobar si esto cumple la última condición, es decir, si a_2\leq\frac{a_1}{1+a_1 P} \ \Rightarrow
% \ P \leq -\frac{1}{a_1}+\frac{1}{a_2}.
% Por tanto, si  P \leq -\frac{1}{a_1}+\frac{1}{a_2}, el punto \mathbf x =(P,0) es KKT.
%
% En resumen:
% 1)   Si  P \geq \frac{1}{a_1}-\frac{1}{a_2},  \ P \geq \frac{1}{a_2}-\frac{1}{a_1}, el punto \mathbf x =\left(\frac{P}{2}-\frac{1}{2 a_1}+\frac{1} {2a_2}, \quad
% \frac{P}{2}+\frac{1}{2 a_1}-\frac{1} {2a_2}\right) es el único punto KKT.
% 2)  Nunca se cumple.
% 3)   Si  P \leq \frac{1}{a_1}-\frac{1}{a_2}, el punto  \mathbf x =(0,P) es el único punto KKT.
% 4)   Si  P \leq -\frac{1}{a_1}+\frac{1}{a_2}, el punto \mathbf x =(P,0) es el único punto KKT.
% Se ve que los casos 1,3,4 son excluyentes, es decir, si se cumple uno no se cumplen los otros dos. Además, una de los casos 1,2,3 se cumple siempre. Por tanto, la solución del problema es el punto KKT correspondiente al caso que se cumple.
%
%% Se analizan algunos casos particulares:
%
% P=1, \ a_1=1 y a_2=2      \Rightarrow \quad \mathbf x^* = \left(\frac{1}{4}\ , \ \frac{3}{4}\right)       (caso 1)
% Se comprueba gráficamente dibujando la función y las restricciones:

a1=1;
a2=2;
P=1;

funcion = @(x1,x2) -log(1+a1*x1) -log(1+a2*x2) ;
restriccion1 = @(x1,x2) x1+x2-P;
restriccion2 = @(x1,x2) -x1;
restriccion3 = @(x1,x2) -x2;
region = [-.5 2 -.5 2];

fcontour(funcion, region, 'LevelStep', 0.1)
colorbar
xlabel('x_1'); ylabel('x_2')
hold on
fimplicit(restriccion1, region,'k')
fimplicit(restriccion2, region,'r')
fimplicit(restriccion3, region,'r')
plot(P/2-1/(2*a1)+1/(2*a2), P/2+1/(2*a1)-1/(2*a2),'ro')
title('f(x), g(x)=0, h1(x)=0, h2(x) = 0')
hold off

%
% P=1, \ a_1=\frac{1}{2} y a_2=2   \Rightarrow \quad \mathbf x^* = (0,1)       (caso 3)
% Se comprueba gráficamente dibujando la función y las restricciones:

a1=1/2;
a2=2;
P=1;

funcion = @(x1,x2) -log(1+a1*x1) -log(1+a2*x2) ;
region = [-.5 2 -.5 2];

fcontour(funcion, region, 'LevelStep', 0.1)
colorbar
xlabel('x_1'); ylabel('x_2')
hold on
fimplicit(restriccion1, region,'k')
fimplicit(restriccion2, region,'r')
fimplicit(restriccion3, region,'r')
plot(0,P,'ro')
title('f(x), g(x)=0, h1(x)=0, h2(x) = 0')
hold off

%
% P=1, \ a_1=1 y a_2=\frac{1}{3}   \Rightarrow \quad \mathbf x^* = (1,0)       (caso 4)
% Se comprueba gráficamente dibujando la función y las restricciones:

a1=1;
a2=1/3;
P=1;

funcion = @(x1,x2) -log(1+a1*x1) -log(1+a2*x2) ;
region = [-.5 2 -.5 2];

fcontour(funcion, region, 'LevelStep', 0.1)
colorbar
xlabel('x_1'); ylabel('x_2')
hold on
fimplicit(restriccion1, region,'k')
fimplicit(restriccion2, region,'r')
fimplicit(restriccion3, region,'r')
plot(P,0,'ro')
title('f(x), g(x)=0, h1(x)=0, h2(x) = 0')
hold off

%
%
%% Métodos numéricos: Función fmincon de Matlab
%
% La función de matlab fmincon resuelve numéricamente problemas de minimización con restricciones, tanto de igualdad como de desigualdad. Utiliza algoritmos iterativos para encontrar los mínimos.
% Las entradas a fmincon son la función a minimizar, f(\mathbf x) ,  las funciones asociadas a las restricciones \left\{h_k(\mathbf x)\right\}, \left\{g_k(\mathbf x)\right\}, y un puntor inicial \mathbf x_0 . A partir del punto inicial el algoritmo empleado comienza a iterar buscando un mínimo.
% Cuando converge a un mínimo a un mínimo, el algoritmo deja de iterar y devuelve el mínimo encontrado.
% La selección del valor inicial suele ser crítica, ya que determina la convergencia del algoritmo y el mínimo encontrado (si es que converge).
% fmincon estima las derivadas de las funciones numéricamente, salvo que se las proporcione explicitamente. En este caso se acelera la convergencia del algoritmo.
% En general, la llamada a la función es:  x = fmincon( f, x0, A, b, Ae, be, lb, ub, fnolcon ),
% donde f es la función objetivo, x0 es el punto semilla, A y b definen las restricciones de desigualdad lineales, Ae y be definen las restricciones de igualdad lineales, lb son las cotas inferiores de las variables,  ub son las cotas superiores de las variables, y fnolcon es una función donde se definen las restricciones no lineales (tanto de igualdad como de desigualdad).
%
%
%% Ejemplo 1:
%                                 \min \quad x_1^2+x_2^2 \\
%  s.t. \quad 1-x_1-x_2 \leq 0 \\
%   \quad \quad 2 x_1-x_2 = 0
% Solución:
% D = \mathbb{R}^2.
% Por tanto,   f(\mathbf x) = x_1^2 + x_2^2, \quad \quad
% h(\mathbf x) = 1 - x_1 - x_2, \quad \quad
% g(\mathbf x) = 2 x_1-x_2
% Las funciones asociadas a la restricciones son lineales. Por tanto, se pueden escribir así:
%  \mathbf A \mathbf x \leq \mathbf b, donde \mathbf A = [- 1,- 1], \ \mathbf b = -1
%
% },
% \mathbf A_e \mathbf x = \mathbf b_e, donde \mathbf A_e = [2,- 1], \ \mathbf b_e = 0
%
% }.
% Busca mínimo:

f = @(x) x(1)^2 + x(2)^2;
A = [-1,-1]; b = -1;
Ae = [2,-1]; be = 0;
x0 = [1,2]; 
x = fmincon(f,x0,A,b,Ae,be,[],[],[]);
disp(x)

%
% Se comprueba gráficamente:

funcion = @(x1,x2) x1.^2 + x2.^2;
restriccion1 = @(x1,x2) 1 - x1 - x2;
restriccion2 = @(x1,x2) 2*x1 - x2;
region = [-1 1 -1 1];

fcontour(funcion, region, 'LevelStep', 0.04)
xlabel('x_1'); ylabel('x_2')
colorbar
hold on
fimplicit(restriccion1, region,'r')
fimplicit(restriccion2, region,'r')
plot(x(1),x(2),'ko')
title('f(x), g(x)=0, h(x)=0')
hold off

%
%
%% Ejemplo 4
%                                 \min \quad -x_1^2 + 2 x_1 + 4 x_2 \\
%  s.t. \quad -x_1 +4 x_2  = 0 \\
% \quad \quad x_1+x_2-4 \leq 0 \\
% \quad \quad -x_1+2 \leq 0
% Solución:
% D = \mathbb{R}^2.

f = @(x) -x(1)^2 + 2*x(1) + 4*x(2);
A = [1,1]; b = 4;
Ae = [-1,4]; be = 0;
lb=[2,-inf];
x0 = [1,2]; 
x = fmincon(f,x0,A,b,Ae,be,lb,[],[]);
disp(x)

%
% Se comprueba gráficamente:

funcion = @(x1,x2) -x1.^2+4*x2+2*x1;
restriccion1 = @(x1,x2) -x1+4*x2;
restriccion2 = @(x1,x2) x2+x1-4;
restriccion3 = @(x1,x2) -x1+2;
region = [0 4 0 4];
fcontour(funcion, region, 'LevelStep', 0.1)
colorbar
xlabel('x_1'); ylabel('x_2')
hold on
fimplicit(restriccion1, region,'r')
fimplicit(restriccion2, region,'r')
fimplicit(restriccion3, region,'r')
plot(x(1),x(2),'ko')
hold off
title('f(x), g(x)=0, h_1(x), h_2(x)')

%
%
%% Ejemplo 5
%                               \min \quad -\log(1+a_1 x_1) - \log(1+a_2 x_2) \\
%  s.t. \quad x_1+x_2 -P= 0 \\
% \quad \quad -x_1 \leq 0 \\
% \quad \quad -x_2 \leq 0
% cuando a_1=1, \ a_2=2, \ P=1.
%
% Solución:
% D = \mathbb{R}^2.

a1=1;
a2=2;
P=1;
f = @(x) -log(1+a1*x(1)) -log(1+a2*x(2));
Ae = [1,1]; be = P;
lb=[0,0];
x0 = [1,1]; 
x = fmincon(f,x0,[],[],Ae,be,lb,[],[]);
disp(x)

%
% Se comprueba gráficamente:

funcion = @(x1,x2) -log(1+a1*x1) -log(1+a2*x2) ;
restriccion1 = @(x1,x2) x1+x2-P;
restriccion2 = @(x1,x2) -x1;
restriccion3 = @(x1,x2) -x2;
region = [-.5 2 -.5 2];
fcontour(funcion, region, 'LevelStep', 0.1)
colorbar
xlabel('x_1'); ylabel('x_2')
hold on
fimplicit(restriccion1, region,'r')
fimplicit(restriccion2, region,'r')
fimplicit(restriccion3, region,'r')
plot(x(1), x(2),'ko')
title('f(x), g(x)=0, h_1(x)=0, h_2(x) = 0')
hold off

%
%
%% Ejemplo 6
%
%                     \min \quad -x_1^3 - x_2^2 - 2 x_1 x_3^2 \\
%  s.t. \quad 2 x_1 + x_2^2 + x_3 - 5 = 0 \\
% \quad \quad -5 x_1^2 + x_2^2 + x_3  \leq -2 \\
% \quad \quad x_1, x_2, x_3 \geq 0
% Solución:
% D = \mathbb{R}^3.

f = @(x) -x(1)^3 - x(2)^2 - 2*x(1)*x(3)^2;
x0 = [1,1,2]; % punto factible
lb = [0,0,0];
x = fmincon(f,x0,[],[],[],[],lb,[],@restricciones_ejemplo6);
disp(x)

%
%
%% Ejemplo 7
%
% Un inversor desea invertir en acciones de dos empresas:  A_1 y  A_2 .
% Las rentablilidades (por euro invertido) de las acciones se modela como variables aleatorias Gaussianas:  Y_i \sim N(\eta_i, \sigma_i^2), \ i=1,2. La varianza \sigma_i^2 se conoce como riesgo de la inversión en A_i y la media \eta_i  es la rentabilidad esperada al invertir en A_i.
% El inversor dispone de P euros para invertir, y debe decidir cuantos euros invierte en cada acción:  x_i, \ i = 1,2.  Por tanto, x_1 + x_2 = P.
% La rentabilidad total vendrá dada por la variable aleatoria  Y = x_1 Y_1 + x_2 Y_2.
% La rentabilidad esperada es la media de Y:    \eta_Y = E[Y] = x_1 \eta_1 + x_2 \eta_2 ,
% El riesgo es la varianza de Y:   \sigma_Y^2 = var [Y] = x_1^2 \sigma_1^2 + x_2^2 \sigma_2^2 + 2 x_1 x_2 \sigma_{1,2} ,
% donde \sigma_{1,2} es la covarianza de las variables aleatorias Y_1 y Y_2. Si las rentabilidades están incorreladas, \sigma_{1,2}=0.
% El inversor busca un compromiso entre la rentabilidad esperada y el riesgo de la inversión.
% Para ello se propone maximizar  la función \alpha \ \eta_Y^2 -  (1-\alpha) \ \sigma_Y^2, donde \alpha  es un coeficiente que refleja la importancia que el inversor da a \eta_Y en relación con \sigma_Y^2.
% Si elige \alpha=1,  el inversor buscará maximizar la rentabilidad esperada sin importarle el riesgo (inversor arriesgado), mientras que cuando elige \alpha = 0, buscará minimizar el riesgo sin importarle la rentabilidad esperada. En la práctica el parámetro toma valores en el intervalo 0 < \alpha < 1
% Para un valor dado de \alpha , y de los momentos \mu_1, \mu_2, \sigma_1^2, \sigma_2^2, \sigma_{1,2}, determine cuantos euros debería invertir en cada acción
%
% Solución:
% D = \mathbb{R}^2.
% La función objetivo será:  f(\mathbf x) = - \alpha \ \eta_Y^2 +  (1-\alpha) \sigma_Y^2 =
% a_1 x_1^2 + a_2 x_2^2 + a_{1,2} x_1 x_2,
% donde a_1 = -\alpha \eta_1^2 + (1-\alpha) \sigma_1^2, \quad
% a_2 = -\alpha \eta_2^2 + (1-\alpha) \sigma_2^2, \quad
% a_{1,2} = -2 \alpha \eta_1 \eta_2 + 2 (1-\alpha) \sigma_{1,2}
% Escribo el problema en la forma estandar
%                     \min \quad a_1 x_1^2 + a_2 x_2^2 + a_{1,2} x_1 x_2  \\
%  s.t. \quad x_1+x_2 -P= 0 \\
% \quad \quad -x_1 \leq 0 \\
% \quad \quad -x_2 \leq 0
%
% Considero el siguiente caso particular:

media1 = 0.2; varianza1 = 0.2;
media2 = 0.3; varianza2 = 0.4;
covarianza = 0.1;
alfa = 0.7;
a1 = -alfa*media1^2+(1-alfa)*varianza1;  
a2 = -alfa*media2^2+(1-alfa)*varianza2;
a12 = -2*alfa*media1*media2 + 2*(1-alfa)*covarianza; 
P = 1000;

%
% Solución con fmincon:

f = @(x) a1*x(1).^2 + a2*x(2).^2 + a12*x(1).*x(2);
Ae = [1,1]; be = P;
lb = [0,0];
x0 = [P/2,P/2]; % punto factible
x = fmincon(f,x0,[],[],Ae,be,lb,[],[]);
disp(x)

%
% Se comprueba gráficamente:

funcion = @(x1,x2) a1*x1.^2 + a2*x2.^2 + a12*x1.*x2;
restriccion1 = @(x1,x2) x1+x2-P;
restriccion2 = @(x1,x2) -x1;
restriccion3 = @(x1,x2) -x2;
region = [-P/10,P+P/10,-P/10,P+P/10];

fcontour(funcion, region, 'LevelStep', 2000)
colorbar
xlabel('x_1'); ylabel('x_2')
hold on
fimplicit(restriccion1, region,'k')
fimplicit(restriccion2, region,'r')
fimplicit(restriccion3, region,'r')
title('f(x), g(x)=0, h_1(x)=0, h_2(x) = 0')
plot(x(1),x(2),'ko')
hold off

%
%
%% Ejercicios propuestos
%
%% Ejercicio 1
%
% Considere el siguiente problema:
%                     \min \quad 1-x_1 +  x_2^2 \\
%  s.t. \quad x_1^2 +  x_2^2 -1 \leq 0 \\
% \quad \quad -2x_1 + x_2 -1 = 0
% Determine los puntos irregulares
% Determine los puntos KKT y los multiplicadores correspondientes
% Determine el punto solución
% Compruebe la solución gráficamente
% Obtenga la solución con fmincon
%
%% Ejercicio 2
%
% Considere el siguiente problema:
%                     \max \quad 4 x_1 + 2 x_2 \\
%  s.t. \quad x_1 - x_2 +1 =0 \\
% \quad \quad x_1 + x_2 \leq 2
% Determine los puntos irregulares
% Determine los puntos KKT y los multiplicadores correspondientes
% Determine el punto solución
% Compruebe la solución gráficamente
% Obtenga la solución con fmincon
%
%% Ejercicio 3
%
% Considere la región de puntos del plano  x_1 + x_2 + x_3 = 1 contenidos dentro de la esfera de radio unidad. Se desea determinar qué punto de la región está más cerca del punto \mathbf p=(1,2,3)  y a qué distancia se encuentra.
% Plantee el problema de optimización
% Resuelvalo con fmincon.
%
%
%% Ejercicio 4
%
% Se desea determinar el cilindro de mínimo volumen con área A = 16, radio r \geq 1 y altura h \geq 1.
% Plantee el problema de optimización
% Determine los puntos irregulares
% Determine los puntos KKT y los multiplicadores correspondientes
% Determine el punto solución
% Compruebe la solución gráficamente
% Obtenga la solución con fmincon
%
%
%% Ejercicio 5
%
% Resuelva el siguiente ejercicio utilizando la función fmincon.
%
%                                         min \quad x_1^2 + 2 x_2^2 + 3 x_3^2  \\
%  s.t. \quad -5 x_1 + x_2 +3 x_3 + 3 \leq 0 \\
% \quad \quad -2 x_1 - x_2 - 2 x_3 +6 = 0  \\
% \quad \quad - x_1^2 + (x_2 + x_3)^2 \leq 1 \\
% \quad \quad x_2 \log x_1 - x_3 \leq 0 \\
% \quad \quad x_1, x_2, x_3 \geq 0 \\
% \quad \quad x_2, x_3 \leq 1
%
%
%% Funciones no lineales de las restricciones:
%

function [c,ceq] = restricciones_ejemplo6(x)
c = -5*x(1).^2 + x(2).^2 + x(3) + 2;
ceq = 2*x(1) + x(2).^2 + x(3) - 5;
end

%