%% Optimización con restricciones de igualdad.
% 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
%                 - Definiciones
%                 - Condiciones necesarias de Lagrange
%                 - Clasificación de los puntos
%                 - Interpretación geométrica
%                 - Método alternativo: sustitución
%                 - Ejemplos
%                 - Ejercicios
%%   - Optimización con varias restricciones de igualdad
%                 - Definiciones
%                 - Condiciones necesarias de Lagrange
%                 - Clasificación de los puntos
%                 - Ejemplos
%                 - Ejercicio
%%   - Métodos numéricos: función fmincon
%                 - Ejemplos
%                 - Ejercicio
%
%
%% Introducción y planteamiento del problema
%
% Determinar los mínimos de una función f(\mathbf x)   (función objetivo) donde la variable independiente \mathbf x está sujeta a restricciones de igualdad de la forma  g_k(\mathbf x) =  0, \ \  k=1,\ldots,K.
% Las funciones f(\mathbf x), \ g_k(\mathbf x) , \ \  k=1,\ldots,K 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_k(\mathbf x) :  D \subset \mathbb{R}^M \longrightarrow \mathbb{R}, \quad k=1,\ldots K.
% Las funciones  f(\mathbf x), \ g_k(\mathbf x) , \ \  k=1,\ldots,K  son derivables. También se supone que K < M.
%
% Notación:
%                         \min \quad f(\mathbf{x}) \\
% \matrix{s.t. \quad g_1(\mathbf x) = 0 \cr \vdots \cr \quad \quad g_K(\mathbf x) = 0}
%
% Los máximos de f(\mathbf{x}) coinciden con los mínimos de  -f(\mathbf{x}).
%
%%
%% Optimización con una restricción de igualdad (K=1)
%
%                 min \quad f(\mathbf{x}) \\
% s.t. \quad g(\mathbf x) = 0
%
%% Definiciones
%
% - Puntos factibles:   Son los puntos que cumplen la restricción,  g(\mathbf{x}) = 0
% - Puntos regulares:  Son los puntos factibles en los que \nabla g(\mathbf x) \ne  \mathbf 0
% .
% - Puntos estacionarios:  Son puntos que cumplen las condiciones de Lagrange:
%
%              \left\{
% \matrix{\nabla f(\mathbf x) +  \lambda  \ \nabla g(\mathbf x) = \mathbf 0 \cr
% g(\mathbf x) = 0}
%
% donde \lambda es un número real que se conoce como multiplicador de Lagrange.
% Las condiciones de Lagrange constituyen un sistema de M+1 ecuaciones con M+1 incógnitas:  \mathbf x, \lambda
%
% A veces las ecuaciones que definen los puntos estacionarios se escriben de forma compacta así:  \nabla_{\mathbf x,\lambda}  L(\mathbf x, \lambda) =  \mathbf 0
% , donde
% L(\mathbf x, \lambda) =  f(\mathbf x) + \lambda \ g(\mathbf x)
%  es la llamada función Lagrangiana.
%
%% Condiciones necesarias de Lagrange
% Si \mathbf x^* es un máximo o mínimo (local o global) y además es un punto regular, entonces es un punto estacionario.
%
%% Clasificación de los puntos
%
% D:
% \left\{
% \matrix{
% - \ \textnormal{Puntos factibles}:
% \left\{\matrix{
% - \ \textnormal{Puntos regulares}:
% \left\{\matrix{
% - \ \textnormal{Puntos estacionarios} ^{(*)} \
% \cr
% - \ \textnormal{Puntos no estacionarios}
% }
% \cr
% \cr
% - \ \textnormal{Puntos irregulares}^{(*)} \quad \quad \quad \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 \quad \quad \quad
% }
%
% (*) Los máximos y los mínimos pertenecen al conjunto de los puntos irregulares o al conjunto de los puntos estacionarios.
% Normalmente, estos conjuntos contienen un número finito de puntos, o ninguno.
%
%% Interpretación geométrica
% Ejemplo: M = 2
%
% Si el punto solución es regular, los vectores gradiente \nabla f(\mathbf x)
%  y \nabla g(\mathbf x)
%  en el punto solución son paralelos.
%
%% Método alternativo: Sustitución
% Partiendo de g(\mathbf x) se intenta despejar una componente de \mathbf x y se sustituye en f(\mathbf x). Con esto el problema se transforma  en un problema de optimización sin restricciones.
% Es menos conveniente que el método de los multiplicadores de Lagrange, salvo en algunos problemas sencillos.
%
%
%% Ejemplo 1
%
% Considere la recta en el plano que pasa por los puntos (0,1) y (1,0). Determine el punto de la recta más próximo al origen.
%
% Solución:
% D = \mathbb{R}^2,
% Ecuación de la recta:   x_2 = -x_1 + 1.  Distancia al origen  \sqrt{x_1^2+x_2^2}
% Por tanto, el problema de optimización es
%                                  \min \quad x_1^2+x_2^2 \\
%  s.t. \quad x_1+x_2-1 = 0
%
% f(\mathbf x) = x_1^2 + x_2^2, \quad \quad g(\mathbf x) = x_1 + x_2 -1
% Si no hubiera restricciones, el mínimo de f(\mathbf x) estaría en  \mathbf x = \mathbf 0. Pero este punto no es factible.
% Los puntos factibles son los de la recta  g(\mathbf x) = x_1 + x_2 -1=0.

% Representa f(x)
funcion = @(x1,x2) x1.^2 + x2.^2;
region = [-1 1 -1 1];
fcontour(funcion, region, 'LevelStep', 0.04)
xlabel('x_1'); ylabel('x_2')
colorbar 

% Representa g(x)
hold on
restriccion = @(x1,x2) x1 + x2 - 1;
fimplicit( restriccion, region,'r')
colorbar
title('f(x), g(x)=0')
hold off

%
% Puntos irregulares:
%         \nabla g(\mathbf x) = (1,1) \ne \mathbf 0
% . Por tanto, no hay puntos irregulares.
%
% Puntos estacionarios:
%          \left\{
% \matrix{\nabla f(\mathbf x) +  \lambda  \ \nabla g(\mathbf x) = \mathbf 0 \cr
% g(\mathbf x) = 0}
% \quad \Rightarrow \quad  \left\{ \matrix{2 x_1 + \lambda = 0 \cr 2 x_2 + \lambda = 0 \cr x_1 + x_2 -1 =0}
% Este es un sistema lineal de ecuaciones,  \mathbf A \mathbf y = \mathbf b, donde
% \mathbf A = \left( \matrix{2 & 0 & 1 \cr 0 & 2 & 1 \cr 1 & 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 1} \right )
% Como \det(\mathbf A) \neq 0, el sistema tiene una única solución:     x_1=\frac{1}{2}, \ \  x_2=\frac{1}{2},  \ \  \lambda = -1.

% Resuelve el sistema de ecuaciones
A=[2 0 1;0 2 1; 1 1 0]; b=[0 0 1]';
y = A\b;

% Dibuja la solucion
fcontour(funcion, region, 'LevelStep', 0.04)
xlabel('x_1'); ylabel('x_2')
colorbar

hold on
fimplicit(restriccion, region,'r')
plot(y(1),y(2),'ro')
hold off
title('f(x), g(x)=0')
disp(['solución: ',num2str(y')])

% Por tanto,  \mathbf x = (1/2, \ 1/2)  es el único punto estacionario. Gráficamente se ve que es el mínimo global.
%
%% Ejemplo 2
%
% Determinar el punto más cercano al origen de entre los puntos de la curva
% x_2^2  = (x_1-1)^3
%
% Solución:
% D = \mathbb{R}^2,
% Distancia al origen  \sqrt{x_1^2+x_2^2}
% Por tanto, el problema de optimización es
%                     \min \quad x_1^2+x_2^2 \\
%  s.t. \quad (x_1-1)^3- x_2^2 = 0
%
% Entonces,    f(\mathbf x) = x_1^2 + x_2^2, \quad \quad
% g(\mathbf x) = (x_1-1)^3- x_2^2
% Representa las funciones anteriores

% Representa f(x)
funcion = @(x1,x2) x1.^2 + x2.^2;
region = [-2 2 -2 2];
fcontour(funcion, region, 'LevelStep', 0.1)
xlabel('x_1'); ylabel('x_2')
colorbar 

% Representa g(x)
hold on
restriccion = @(x1,x2) (x1-1).^3 - x2.^2;
fimplicit( restriccion, region,'r')
colorbar
title('f(x), g(x)=0')
hold off

% Puntos irregulares:
%  \nabla g(\mathbf x) =\left[ \matrix{2 (x_1-1)^2 \cr -2 x_2} \right]. Por tanto el punto (1, 0) es un punto irregular (ver que es factible).  La gráfica muestra que es un mínimo.
% Puntos estacionarios:
%             \left\{
% \matrix{\nabla f(\mathbf x) +  \lambda  \ \nabla g(\mathbf x) = \mathbf 0 \cr
% g(\mathbf x) = 0}
% \quad \Rightarrow \quad  \left\{ \matrix{
% 2 x_1 + 2 \lambda (x_1-1)^2 = 0 \cr
% 2 x_2 - 2 \lambda x_2 = 0 \cr
% (x_1-1)^3 - x_2^2 = 0}
% Es fácil comprobar que este sistema no tiene solución. Por tanto, no hay puntos estacionarios.
% Este es un ejemplo en el que la solución del problema es un punto irregular.
%
%% Ejemplo 3
%
%                     \min \quad x_1^3+x_2^3 \\
%  s.t. \quad x_1- x_2 = 0
%
% En este ejemplo  D = \mathbb{R}^2,    f(\mathbf x) = x_1^3 + x_2^3, \quad \quad
% g(\mathbf x) = x_1 - x_2

% Representa f(x)
funcion = @(x1,x2) x1.^3 + x2.^3;
region = [-2 2 -2 2];
fcontour(funcion, region, 'LevelStep', 0.1)
xlabel('x_1'); ylabel('x_2')
colorbar 

% Representa g(x)
hold on
restriccion = @(x1,x2) x1 - x2;
fimplicit( restriccion, region,'r')
colorbar
title('f(x), g(x)=0')
hold off

% Gráficamente se ve que no hay máximos ni mínimos.
%
% Puntos irregulares:
%  \nabla g(\mathbf x) =\left[ \matrix{
% 1 \cr
% -1
% } \right]. Por tanto no hay puntos irregulares.
%
% Puntos estacionarios:
%   \left\{
% \matrix{\nabla f(\mathbf x) +  \lambda  \ \nabla g(\mathbf x) = \mathbf 0 \cr
% g(\mathbf x) = 0}
% \quad \Rightarrow \quad  \left\{ \matrix{
% 3 x_1^2 + \lambda = 0 \cr
% 3 x_2^2 - \lambda = 0 \cr
% x_1 - x_2 = 0}
% Se resuelve a mano facilmente. Sólo hay un punto estacionario:  \mathbf x = ( 0, 0). Este punto no es máximo ni mínimo.
% Este ejemplo muestra que las condiciones de Lagrange no son suficientes. Puede haber puntos estacionarios que no sean máximos ni mínimos.
%
%% Ejemplo 4
%
% De entre los triangulos rectángulos de area A, determine el que tiene menor hipotenusa.
%
%  Solución:
% x_1: primer cateto,    x_2:segundo cateto
% D =  \{ \mathbf x = (x_1,x_2) | x_1,x_2 > 0 \} \subset \mathbb{R}^2
% Hipotenusa: \sqrt{x_1^2 + x_2^2}
% Area del triángulo:  A = \frac{x_1 \ x_2}{2}
% Problema de optimización:
%                     \min \quad x_1^2 +  x_2^2 \\
%  s.t. \quad x_1  x_2 -2 A=0
%
% Por tanto,   f(\mathbf x) =  x_1^2 + x_2^2, \quad  g(\mathbf x)= x_1 x_2 - 2A
%
% Puntos irregulares:
%         \nabla g(\mathbf x) =\mathbf 0  solo se cumple en \mathbf x = \mathbf 0 \notin D.  Por tanto no hay puntos irregulares.
%
% Puntos estacionarios:
%          \left\{
% \matrix{\nabla f(\mathbf x) +  \lambda  \ \nabla g(\mathbf x) = \mathbf 0 \cr
% g(\mathbf x) = 0}
% \quad \Rightarrow \quad  \left\{ \matrix{2 x_1 + \lambda x_2= 0 \cr 2 x_2 + \lambda x_1 = 0 \cr x_1 x_2 - 2 A =0}
% Sistema de ecuaciones no lineal, pero se puede resolver analiticamente. Sólo hay un punto estacionario:
% x_1 = x_2 = \sqrt{2 A}, \  \lambda = -2
% Se comprueba gráficamente para un valor concreto de A:

% Representa f(x) para un valor concreto de A
A=2;
funcion = @(x1,x2) x1.^2 + x2.^2;
region = [-4 4 -4 4];
fcontour(funcion, region, 'LevelStep', 0.4)
xlabel('x_1'); ylabel('x_2')
colorbar
% Representa g(x)
hold on
fimplicit(@(x1,x2) x1.*x2-2*A , region,'r')
title('f(x), g(x)=0')
hold off

% Resuelve el sistema de ecuaciones (no lineal) de los puntos estacionarios numéricamente.

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))])

% Se comprueba gráficamente la solución obtenida.
%
%% Ejemplo 5
%
% Considere la esfera de radio R centrada en el origen, y sea \mathbf p = (p_1,p_2,p_3) un punto en el espacio distinto del origen. Determine el punto de la esfera más cercano a \mathbf p.
%
%  Solución:
% Puntos de la esfera:   x_1^2+x_2^2+x_3^2 = R^2
% Distancia al punto \mathbf p:  \sqrt{(x_1-p_1)^2+(x_2-p_2)^2+(x_3-p_3)^2}
% D =  \mathbb{R}^3
% El problema de optimización:
%                     \min \quad  (x_1-p_1)^2+(x_2-p_2)^2 + (x_3-p_3)^2 \\
%  s.t. \quad x_1^2+x_2^2+x_3^2 = R^2
%
% Por tanto,   f(\mathbf x) = (x_1-p_1)^2+(x_2-p_2)^2+(x_3-p_3)^2 , \quad \quad
% g(\mathbf x) = x_1^2+x_2^2+x_3^2 - R^2
%
% Puntos irregulares:
%  \nabla g(\mathbf x) =\left[ \matrix{
% 2 x_1 \cr
% 2 x_2 \cr
% 2 x_3
% } \right]. Por tanto no hay puntos irregulares (R>0).
%
% Puntos estacionarios:
%   \left\{
% \matrix{\nabla f(\mathbf x) +  \lambda  \ \nabla g(\mathbf x) = \mathbf 0 \cr
% g(\mathbf x) = 0}
% \quad \Rightarrow \quad  \left\{ \matrix{
% 2(x_1-p_1) + \lambda 2 x_1 = 0 \cr
% 2(x_2-p_2) + \lambda 2 x_2 = 0 \cr
% 2(x_3-p_3) + \lambda 2 x_3 = 0 \cr
% x_1^2 + x_2^2 + x_3^2 = R^2}
% Ver que \lambda \ne -1, ya que \lambda = -1 \ \Rightarrow \ p_1=p_2=p_3=0. Entonces el sistema se puede escribir así
%       \left\{ \matrix{
% x_1 = \frac{p_1}{1+\lambda} \cr
% x_2 = \frac{p_2}{1+\lambda} \cr
% x_3 = \frac{p_3}{1+\lambda} \cr
% x_1^2 + x_2^2 + x_3^2 = R^2}
% Sustituyendo las tres primeras ecuaciones en la cuarta, se tiene,  1+\lambda = \pm \frac{\sqrt{p_1^2+p_2^2+p_3^2}}{R}, que si se sustituye en las primeras ecuaciones se obtienen estos dos puntos estacionarios:
% x_1 = \frac{p_1 R}{\sqrt{p_1^2+p_2^2+p_3^2}}, \quad
% x_2 = \frac{p_2 R}{\sqrt{p_1^2+p_2^2+p_3^2}}, \quad
% x_3 = \frac{p_3 R}{\sqrt{p_1^2+p_2^2+p_3^2}}, \quad
% \lambda = -1 + \frac{\sqrt{p_1^2+p_2^2+p_3^2}}{R}
% x_1 = \frac{-p_1 R}{\sqrt{p_1^2+p_2^2+p_3^2}}, \quad
% x_2 = \frac{-p_2 R}{\sqrt{p_1^2+p_2^2+p_3^2}}, \quad
% x_3 = \frac{-p_3 R}{\sqrt{p_1^2+p_2^2+p_3^2}}, \quad
% \lambda = -1 - \frac{\sqrt{p_1^2+p_2^2+p_3^2}}{R}
% Uno de los puntos estacionarios es un mínimo (punto de la esfera más cercano a \mathbf p) y el otro es un máximo (punto de la esfera más alejado de \mathbf p)
%
% A continuación se comprueba un caso particular resolviendo el sistema de ecuaciones (no lineal) numéricamente.

R=1;
p=[1,1.5,-1];
sistema = @(y)[2*(y(1)-p(1))+2*y(1)*y(4) , 2*(y(2)-p(2))+2*y(2)*y(4) , 2*(y(3)-p(3))+2*y(3)*y(4) ,y(1)^2+y(2)^2+y(3)^2-R^2];
y0=[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))])

% Ver que fsolve converge a una de las soluciones en función del puntos semilla elegido y0.
%
%% Ejercicio 1
%
% Aplique el método de los multiplicadores de Lagrange para obtener el punto que minimiza de la función  f(\mathbf x) = (x_1-2)^2 + (x_2-1)^2 + (x_1-2) (x_2-1)
% de entre los puntos que pertenecen a la recta de pendiente unidad que pasa por el origen.
% Plantee el problema como un problema de optimización
% Determine los puntos irregulares
% Plantee el sistema de ecuaciones que cumplen los puntos estacionarios
% Obtenga los puntos estacionarios manualmente y numéricamente.
%
%% Ejercicio 2
%
% Sea el elipsoide  x_1^2 + 4  x_2^2 + 9 x_3^2 = 1.  Se desea determinar el punto del elipsoide más cercano al punto  \mathbf p =(1, -1, 2).
% Plantee el problema como un problema de optimización
% Determine los puntos irregulares
% Plantee el sistema de ecuaciones que cumplen los puntos estacionarios
% Determine los puntos estacionarios.
% Determine el punto solución
%
%% Ejercicio 3
%
% Se pretende determinar el cilindro de menor superficie de entre todos los cilindros de volumen V=10.
% Plantee el problema como un problema de optimización
% Resuelva el problema con el método de sustitución
% Determine los puntos irregulares
% Determine los puntos estacionarios
% Determine la solución del problema y compruebela gráficamentea aplicando el método de los multiplicadores de Lagrange y compruebe el resultado gráficamente.
%
%%
%% Optimización con varias restricciones de igualdad
%
%                 \min \quad f(\mathbf{x}) \\
% s.t. \quad g_k(\mathbf x) = 0, \quad k=1, \ldots, K
%
%% Definiciones:
%
% - Puntos factibles:   Son los puntos que cumplen las restricciones,  g_k(\mathbf{x}) = 0, \quad k=1, \ldots, K.
% - Puntos regulares:  Son los puntos factibles en los que los vectores gradiente \{\nabla g_k(\mathbf x)\}_{k=1}^K son linealmente independientes.
% - Puntos estacionarios:  Son puntos regulares que cumplen las condiciones de Lagrange:
%
%
% \left\{ \matrix {\nabla f(\mathbf x) + \sum_{k=1}^K \lambda_k \ \nabla g_k(\mathbf x) = \mathbf 0 \cr
%   g_1(\mathbf x) = 0 \cr
% \vdots \cr
% g_K(\mathbf x) = 0
% }
%
%
% donde \lambda_1, \ldots, \lambda_K son números reales llamados multiplicadores de Lagrange.
% Las condiciones de Lagrange constituyen un sistema de M+K ecuaciones con M+K incógnitas:  \mathbf x, \lambda_1, \ldots, \lambda_K
%
%% Condiciones necesarias de Lagrange:
% Si \mathbf x^* es un máximo o mínimo (local o global) y además es un punto regular, entonces es un punto estacionario.
%
%% Clasificación de los puntos:
%
% D:
% \left\{
% \matrix{
% - \ \textnormal{Puntos factibles}:
% \left\{\matrix{
% - \ \textnormal{Puntos regulares}:
% \left\{\matrix{
% - \ \textnormal{Puntos estacionarios} ^{(*)} \
% \cr
% - \ \textnormal{Puntos no estacionarios}
% }
% \cr
% \cr
% - \ \textnormal{Puntos irregulares}^{(*)} \quad \quad \quad \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 \quad \quad \quad
% }
%
% (*) Los máximos y los mínimos pertenecen al conjunto de los puntos irregulares o al conjunto de los puntos estacionarios.
% Normalmente, estos conjuntos contienen un número finito de puntos, o ninguno.
%
%
%% Ejemplo 6
%
%                              \min \quad  -x_1^2  - x_2^2  - 2 x_1 x_2 + x_1 x_3 + x_2 x_4  \\
%  \quad s.t. \quad 2 x_1 + x_3 = 12 \\
% \quad \quad \quad 4 x_2 + x_4 = 32
%
% Solución:
% En este ejemplo   f(\mathbf x) = -x_1^2  - x_2^2  - 2 x_1 x_2 + x_1 x_3 + x_2 x_4, \quad \quad g_1(\mathbf x) = 2 x_1 + x_3 - 12, \quad g_2(\mathbf x) = 4 x_2 + x_4 - 32
%   D = \mathbb{R}^4.
% Los puntos factibles son los puntos que cumplen:  g_1(\mathbf x) = 2 x_1 + x_3 - 12 = 0, \quad
% g_2(\mathbf x) = 4 x_2 + x_4 - 32 = 0.
% Todos los puntos son puntos regulares porque \nabla g_1(\mathbf x) = [2, 0, 1, 0], \quad  \nabla g_2(\mathbf x) = [0, 4, 0, 1], que son linealmente independientes:

dg1 = [2, 0, 1, 0]'; dg2 = [0, 4, 0, 1]';
disp(rank([dg1 dg2]))

% Por tanto, si hay puntos solución, estos son puntos estacionarios:
%
% \left\{ \matrix {\nabla f(\mathbf x) + \lambda_1 \nabla g_1(\mathbf x) + \lambda_2 \nabla g_2(\mathbf x) = \mathbf 0 \cr
%   g_1(\mathbf x) = 0 \cr
%
% g_2(\mathbf x) = 0
% }
%  \quad \Rightarrow \quad  \left\{ \matrix{ -2 x_1 - 2 x_2 + x_3 + 2 \lambda_1 = 0 \cr
%                       -2 x_1 - 2 x_2 + x_4 + 4 \lambda_2  = 0 \cr
%                      x_1 +  \lambda_1 = 0 \cr
%               x_2 + \lambda_2 =0 \cr
% 2 x_1 + x_3 -12 \cr
% 4 x_2 + x_4 -32}
% Este es un sistema de ecuaciones lineales con una solución:

A=[-2 -2  1  0  2  0 ;...
   -2 -2  0  1  0  4 ;...
    1  0  0  0  1  0 ;...
    0  1  0  0  0  1 ;...
    2  0  1  0  0  0 ;...
    0  4  0  1  0  0]; 
disp(det(A))

% Calculo el punto estacionario:

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

%
% El punto estacionario es un máximo.
%
%
%% Ejemplo 7
%
%                     \min \quad x_3 \\
%  \quad s.t. \quad x_1^2+x_2^2 = 4 \\
% \quad \quad \quad x_1 + x_2 + x_3  = 5
% Solución:
% Los puntos factibles son los puntos que cumplen:  g_1(\mathbf x) = x_1^2 + x_2^2 - 4 = 0, \quad
% g_2(\mathbf x) = x_1 + x_2 + x_3 - 5 = 0.   D = \mathbb{R}^3.
% Puntos regulares:
% \nabla g_1(\mathbf x) = [2 x_1, 2 x_2, 0], \quad  \nabla g_2(\mathbf x) = [1,1,1]. En todos los puntos los vectores gradiente son linealmente independientes salvo \mathbf x = \mathbf 0. Pero este no es un punto factible. Por tanto no hay puntos irregulares. Si hay puntos solución, estos deben ser puntos estacionarios:
%
% \left\{ \matrix {\nabla f(\mathbf x) + \lambda_1 \nabla g_1(\mathbf x) + \lambda_2 \nabla g_2(\mathbf x) = \mathbf 0 \cr
%   g_1(\mathbf x) = 0 \cr
%
% g_2(\mathbf x) = 0
% }
%  \quad \Rightarrow \quad  \left\{ \matrix{ 2 x_1  \lambda_1 + \lambda_2= 0 \cr
%                       2 x_2  \lambda_1 + \lambda_2= 0 \cr
%                      1 + \lambda_2 = 0 \cr
%            x_1^2 + x_2^2 - 4 =0 \cr
%  x_1 + x_2 + x_3-5 =0}
% que es un sistema de ecuaciones no lineales. Calcula numericamente los puntos estacionarios:

funcion = @(x)[2*x(1)*x(4)+x(5), 2*x(2)*x(4)+x(5), 1+x(5), x(1)^2+x(2)^2-4, x(1)+x(2)+x(3)-5];
x0=ones(1,5);
x = fsolve(funcion,x0);
disp(['solución: ',num2str(x)])
disp(['valor de la función en la solución: ',num2str(funcion(x))])

% El sistema de ecuaciones se puede resolver analíticamente y se obtienen dos puntos estacionarios:  \left (\sqrt 2, \sqrt 2, 5-2\sqrt 2,\frac{1}{2 \sqrt 2}, -1 \right), \quad
% \left (-\sqrt 2, -\sqrt 2, 5+\sqrt 2,\frac{-1}{2 \sqrt 2}, -1 \right).
% Por tanto, teniendo en cuenta que f(\mathbf x) = x_3,  el segundo punto estacionario es un máximo y el primero es un mínimo .
%
%
%% Ejercicio 4
%
%                     \min \quad x_1^2 + x_2^2 + x_3^2  \\
%  \quad s.t. \quad x_1 + x_2 = 1 \\
% \quad \quad \quad x_1-x_3 = 0
% Resuelva el problema con el método de sustitución
% Resuelva el problema con el método de los multiplicadores de Lagrange determinando los puntos irregulares y los puntos estacionarios.
%
%
%% Métodos numéricos: Función fmincon de Matlab
%
% La función de matlab fmincon resuelve numéricamente problemas de minimización con restricciones de igualdad. 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 de igualdad g_k(\mathbf x) , 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, 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.
% Trata de forma diferente las restricciones asociadas a funciones g_k(\mathbf x)  lineales y no lineales.
%
%
%% Ejemplo 8:
%
%                     min \quad  x_1^2  + x_2^2 + 2 x_1 x_2 - x_1 x_3 - x_2 x_4  \\
%  \quad s.t. \quad 2 x_1 + x_3 = 12 \\
% \quad \quad \quad 4 x_2 + x_4 = 32
% En este ejemplo   f(\mathbf x) = x_1^2  + x_2^2 + 2 x_1 x_2 - x_1 x_3 - x_2 x_4, \quad \quad g_1(\mathbf x) = 2 x_1 + x_3 - 12, \quad g_2(\mathbf x) = 4 x_2 + x_4 - 32,  D = \mathbb{R}^4.
% Las funciones asociadas a las restricciones son lineales. Por tanto, se puede escribir así: \mathbf A \mathbf x = \mathbf b, donde \mathbf A = \left(\matrix{
%  2 &0 &1&0 \cr
% 0& 4 &0 &1
% } \right), \
% \ \mathbf b = \left(\matrix{
% 12 \cr 32
% } \right).

funcion = @(x) x(1)^2 + x(2)^2 + 2*x(1)*x(2) - x(1)*x(3) - x(2)*x(4) ;
A = [2,0,1,0; 0,4,0,1]; b = [12; 32];
x0 = [3,4,6,16]; % punto factible
x = fmincon(funcion,x0,[],[],A,b,[],[],[]);
disp(x)

%
%% Ejemplo 9
%
%                     \min \quad x_1^2 +  x_2^2 \\
%  s.t. \quad x_1  x_2 -4=0
% Es el mismo problema que en el ejemplo 4.
%
% Solución:
% En este ejemplo   f(\mathbf x) =  x_1^2 + x_2^2, \quad  g(\mathbf x)= x_1 x_2 - 4,   D = \mathbb{R}^2.
% La restricción es no lineal. La función asociada se define en una función Matlab  (restriccion1_ejemplo9) al final del LiveScript. Las funciones Matlab deben ir final del LiveScript.

funcion = @(x) x(1)^2 + x(2)^2 ;
x0 = [4,1]; 
x = fmincon(funcion,x0,[],[],[],[],[],[],@restriccion_ejemplo9);
disp(x)

%
%%
%% Ejemplo 10:
%                     \min \quad x_3 \\
%  \quad s.t. \quad x_1^2+x_2^2 = 4 \\
% \quad \quad \quad x_1 + x_2 + x_3  = 5
%  g_1(\mathbf x) = x_1^2 + x_2^2 - 4 = 0, \quad
% g_2(\mathbf x) = x_1 + x_2 + x_3 - 5 = 0.
% La segunda restricción es lineal:  \mathbf A \mathbf x = \mathbf b, donde \mathbf A = [ 1, 1, 1], \ \mathbf b = 5
%
% }.
% La primera restricción es no lineal. La función asociada se define en una función Matlab. Las funciones Matlab deben ir al final del LiveScript:

funcion = @(x) x(3) ;
A = [1,1,1]; b = 5;
x0 = [0,2,3]; % punto factible
x = fmincon(funcion,x0,[],[],A,b,[],[],@restriccion1_ejemplo10);
disp(x)

% La restricción no lineal  (restriccion1_ejemplo10) está definida al final del livescript
%
%% Ejercicio 5
%
% Se desea determinar el punto de la circunferencia unidad más cercano al punto \mathbf p =(1,-1.5)
% Plantee el problema como un problema de optimización
% Resuelva el problema con la función fmincon.
%

function [c,ceq] = restriccion1_ejemplo10(x)
c = [];
ceq = x(1)^2 + x(2)^2 - 4;
end

%

function [c,ceq] = restriccion_ejemplo9(x)
c = [];
ceq = x(1)*x(2) - 4;
end

%