% Activate the following line by filling in the right side. If for example the name of the root file is Main.tex, write
% "...root = Main.tex" if the chapter file is in the same directory, and "...root = ../Main.tex" if the chapter is in a subdirectory.
 
%!TEX root = ../Metodos_Numericos.tex  

\chapter[Optimización]{Optimización}


\section{Optimización no-restringida. Método de búsqueda de la sección dorada}

La búsqueda de la sección dorada es una técnica simple de búsqueda de una sola variable de propósito general. La clave para hacer eficiente este procedimiento es la mejor elección de los puntos intermedios. Esta meta se puede alcanzar al especificar que las siguientes dos condiciones se cumplan (ver \cite{Chapra00}).

\begin{equation}
l_0 = l_1 + l_2 \label{eq:dor1}
\end{equation}

\begin{equation}
\frac{l_1}{l_0} =\frac{l_2}{l_1}\label{eq:dor2}
\end{equation}

Sustituyendo la ecuación \ref{eq:dor1} en la ecuación
\ref{eq:dor2} obtenemos:
\begin{equation}
\frac{l_1}{l_1 + l_2} =\frac{l_2}{l_1}\label{eq:dor3}
\end{equation}

Definimos $R=\frac{l_2}{l_1}$ y sustituimos en la ecuación \ref{eq:dor3} para obtener $R^2+R-1$, cuya solución positiva es $R=\frac{\sqrt{5} - 1}{2}$. $R$ es definido como la razón dorada y fue un número ampliamente utilizado por los Griegos en su arquitectura.

El algoritmo de la razón dorada es:

\begin{algorithm}
  \caption{Algoritmo de Búsqueda con Razón Dorada}
  \label{alg:rdorada}
  \begin{algorithmic}
  \REQUIRE $x_u, x_l, f(x)$
  \ENSURE $x$
   \WHILE {$x_u - x_l > Tol$}
	\STATE $d =  R(x_u-x_l) $
	\STATE $x_1 = x_l + d $
	\STATE $x_2 = x_u - d $
	\STATE $k = k+1 $
	\IF {$f(x_1) < f(x_2)$}
    		\STATE $x_l = x_2$
    	\ELSE 
		\STATE $x_u = x_1$
    	\ENDIF
   \ENDWHILE
   \IF {$f(x_1) < f(x_2)$} 
   \RETURN $x_1$
   \ELSE 
   	\RETURN {$x_2$}
   \ENDIF
  \end{algorithmic}
\end{algorithm}

La implementación en Matlab es:

\begin{verbatim}
function y = Razon_Dorada(f, xl, xu)

R =  (sqrt(5) -1)/2
t = 0;

while (xu-xl > 0.0001)
    d = R*(xu - xl);
    x1 = xl + d;
    x2 = xu - d;
    t = t+1;
    if f(x1) < f(x2)
        xl = x2;
    else 
        xu = x1;
    end;
    
    if f(x1) < f(x2)
        y = x1;
    else
        y=x2;
   end;    
end;
\end{verbatim}

\subsection{Ejemplo 1}
Use la búsqueda de la sección dorada para encontrar el mínimo de la función $f(x)=-2seno(x)+\frac{x^2}{10}$ en el intervalo [0,4].

Los resultados de este ejemplo se muestran en la siguiente tabla 
  \begin{center}
     \begin{tabular}[b]{|r|r|r|r|r|r|r|r|r|r|}
     \hline
     \textbf{$k$}&\textbf{$x_l$}&\textbf{$f_l$}&\textbf{$x_2$}&\textbf{$f_2$}&\textbf{$x_1$}&\textbf{$f_1$}&\textbf{$x_u$}&\textbf{$f_u$}&\textbf{$d$}\\
     \hline
     0 & 0.0000 & 0.0000 & 1.5279 & -1.7647 & 2.4721 & -0.6300 & 4.0000 & -3.1136 & 2.4721 \\
     1 & 0.0000 & 0.0000 & 0.9443 & -1.5310 & 1.5279 & -1.7647 & 2.4721 & -0.6300 & 1.5279\\
     2 & 0.9443 & -1.5310 & 1.5279 & -1.7647 & 1.8885 & -1.5432 & 2.4721 & -0.6300 &
     0.9443\\
     3 & 0.9443 & -1.5310 & 1.3050 & -1.7595 & 1.5279 & -1.7647 & 1.8885 & -1.5432 &
     0.5836\\
     4 & 1.3050 & -1.7595 & 1.5279 & -1.7647 & 1.6656 & -1.7136 & 1.8885 & -1.5432 &
     0.3607\\
     ... & ... & ... & ... & ... & ... & ... & ... & ... & ... \\
     17 & 1.4267 & -1.7757 & 1.4271 & -1.7757 & 1.4274 & -1.7757 & 1.4278 & -1.7757 & 0.0007\\
     \hline
    \end{tabular}
   \end{center}

Según los resultados de la tabla la solución es $x=1.4267$ con $f(x)=-1.7757$.

Para realizar la ejecución dar

\begin{verbatim}
>> Razon_Dorada(@fun_1, 0, 4)

ans =

    1.4276
\end{verbatim}


\subsection{Ejemplo 2}
Use la búsqueda de la sección dorada para encontrar el mínimo de la función $f(x)=-x e^{-x^2}$ en el intervalo [0,4].

Los resultados de este ejemplo se muestran en la siguiente tabla

  \begin{center}
     \begin{tabular}[b]{|r|r|r|r|r|r|r|r|r|}
     \hline
     \textbf{$k$}&\textbf{$x_l$}&\textbf{$f_l$}&\textbf{$x_2$}&\textbf{$f_2$}&\textbf{$x_1$}&\textbf{$f_1$}&\textbf{$x_u$}&\textbf{$f_u$}\\
     \hline
         0     &   0.0000     &    0.0000  &  1.5279  & -0.1480  &  2.4721 &  -0.0055  &  4.0000 &  -0.0000\\
         1    &    0.0000     &    0.0000  &  0.9443  & -0.3871  &  1.5279 &  -0.1480  &  2.4721 &  -0.0055\\
         2    &    0.0000     &    0.0000  &  0.5836  & -0.4151  &  0.9443 &  -0.3871  &  1.5279  & -0.1480\\
         3    &    0.0000     &    0.0000  &  0.3607  &  -0.3167 &   0.5836 &  -0.4151 &   0.9443 &  -0.3871\\
        4     &    0.3607     &   -0.3167  &  0.5836  &  -0.4151 &   0.7214 &  -0.4287 &   0.9443 &  -0.3871\\
        ... & ... & ... & ... & ... & ... & ... & ... & ...  \\
        22  &  0.7071 &  -0.4289  &  0.7071 &  -0.4289  &  0.7071 &  -0.4289  &  0.7072  & -0.4289 \\
        \hline
   \end{tabular}
   \end{center}

Según los resultados de la tabla la solución es $x=0.7071$ con $f(x)=-0.4289$.

Para realizar la ejecución dar

\begin{verbatim}
>> Razon_Dorada(@fun_2, 0, 4)

ans =

    0.7071
\end{verbatim}

\section{Optimización no-restringida. Método de Newton}

\subsection{Método de Newton en una dimensión}
Consideremos que la función $f: \mathbb{R} \rightarrow \mathbb{R}$ es de clase dos, es decir que la segunda derivada puede ser calculada. La idea  consiste en reemplazar en la vecindad del punto $x^k$ de la función $f$ por una aproximación cuadrática $q(x)$ dada por
$$
q(x)=f(x^{(k)})+  f'(x^{(k)})(x-x^{(k)})+ \frac{1}{2} f''(x^{(k)})(x-x^{(k)})^2
$$
llamaremos a $x^{(k+1)}$ el mínimo de $q(x)$. Para calcularlo, es necesario que la matriz $ f''(x^{(k)})$ sea positivo. La función $q(x)$ es entonces estrictamente convexa y tiene un mínimo único en $x_{k+1}$ dado por
$$
q'(x^{(k+1)})=0
$$
Esto da lugar sistema lineal de ecuaciones:
$$
f'(x^{(k)})=-f''(x^{(k)})(x^{(k+1)}-x^{(k)})
$$
La solución para $x^{(k+1)}$ es el mínimo, lo cual da lugar a la
recurrencia
\begin{equation} \label {eq:newton}
x^{(k+1)}=x^{(k)}-\frac{f'(x^{(k)})}{ f''(x^{(k)})}
\end{equation}

\begin{algorithm}
  \caption{Algoritmo de Newton Unidimensional}
  \label{alg:Newton}
  \begin{algorithmic}
  \REQUIRE $ f(x), x_0$
  \ENSURE $x^{(k)}$
    \WHILE {$\|f'(x^{(k)}) \| \geq \varepsilon$}
     \STATE $x^{(k+1)}=x^{(k)}-\frac{f'(x^{(k)})}{ f''(x^{(k)})} $
     \STATE $k=k+1$
    \ENDWHILE    
    \RETURN $x^{(k)} $
  \end{algorithmic}
\end{algorithm}


La implementación en Matlab es

\begin{verbatim}
function y = Newton_1d(derf, der2f, x)
x1 = x;
x2 = x;
while 1
    delta = derf(x1)/der2f(x1);
    x2 = x1 - delta;
    disp([x1, derf(x1), der2f(x1), x2]);
    if abs(delta) < 1e-6 break;
    else
        x1 =x2;
    end;
end;
y = x2;
\end{verbatim}

\subsection{Ejemplo 1}
Calcular el máximo de la función $f(x)=-2seno(x)+\frac{x^2}{10}$ utilizando el método de Newton con un valor inicial $x_0 = 0.5$.

Para resolver tenemos que calcular $f'(x)$ y $f''(x)$
\[
f'(x) =  \frac{x}{5}  - 2 cos(x) 
\]

\[
f''(x) = \frac{1}{5}+ 2 sin(x) 
\]

Para resolver hacemos

\begin{eqnarray*}
	k &=& 0 \\
 	x^{(1)} & = & 0.5000 - \frac{-1.6551}{1.1588}  = 1.9282  \\
	k &=& 1 \\
	x^{(2)} & = & 1.9282 -\frac{1.0854}{2.0735} = 1.4047 \\
	k &=& 2 \\
	x^{(3)} & = & 1.4047 - \frac{-0.0495}{2.1725} = 1.4275 \\
	k &=& 3 \\
	x^{(4)} & = & 1.4275 - \frac{ 8.20126 \times 10^{-5}}{2.1795}  =1.4275 \\
	k &=& 4 \\
	x^{(5)} & = & 1.4275 -\frac{2.02094\times 10^{-10}}{2.1795} =1.4275\\
\end{eqnarray*}


Para ejecutar hacer 

\begin{verbatim}
>> Newton_1d(@der_fun_1, @der2_fun_1, 0.5)

ans =

    1.4276
    
function y = der_fun_1(x)
     y = x/5 - 2*cos(x);    
     
function y = d2f1(x)
     y = 1/5 + 2*sin(x);     
     
\end{verbatim}

\subsection{Ejemplo 2}

Dada la función $f(x)=-x e^{-x^2}$, calcular el mínimo utilizando el método de Newton para un valor inicial $x_0 = -1.0$.

Para resolver tenemos que calcular $f'(x)$ y $f''(x)$
\[
f'(x) =  (2 x^2 - 1 )e^{-x^2}
\]

\[
f''(x) = (6 x - 4 x^3) e^{-x^2} 
\]
Para soluciona hacemos 
\begin{eqnarray*}
   k &=& 0 \\
   x^{(1)} &=&  -1.0000  - \frac{0.3679}  {-0.7358} = -0.5000\\
   k &=& 1 \\
   x^{(2)} &=&  -0.5000 - \frac{ -0.3894} {-1.9470} =  -0.7000\\
   k &=& 2 \\
   x^{(3)} &=&  -0.7000 - \frac{-0.0123}{-1.7325} =  -0.7071\\
   k &=& 3 \\
   x^{(4)} &=&  -0.7071 - \frac{-0.0001}{-1.7156} =  -0.7071\\  
   k &=& 4 \\
   x^{(5)} &=&  -0.7071 - \frac{-0.0000}{-1.7155} =  -0.7071
\end{eqnarray*}

 
Para ejecutar en Matlab 

\begin{verbatim}
>> Newton_1d(@der_fun_2, @der2_fun_2, -1.0)

ans =

    -0.7071
    
function y = der_fun_2(x)
    y = (2*x^2 -1)*exp(-x^2);

function y = der2_fun_2(x)
    y = (6*x-4*x^3)*exp(-x^2);	
\end{verbatim}


\subsection{Método de Newton en $N$ dimensiones}

Consideremos que la función $f$ es de clase dos, es decir que la segunda derivada puede ser calculada. La idea  consiste en reemplazar en la vecindad del punto $x^k$ de la función $f$ por una aproximación cuadrática $q(x)$ dada por
$$
q(x)=f(x^{(k)})+ \nabla f^T(x^{(k)})(x-x^{(k)})+ \frac{1}{2}(x-x^{(k)})^T \nabla^2
f(x^{(k)})(x-x^{(k)})
$$
donde $\nabla f(x^{(k)})$ es el gradiente de la función o vector de primeras derivadas y $\nabla f^2(x^{(k)})$ es la matriz de segundas derivadas o matriz Hessiana, ambos valuados en el muno $x^{(k)}$. Llamaremos a $x^{(k+1)}$ el mínimo de $q(x)$ y para calcularlo, es necesario que la matriz $\nabla^2 f(x^{(k)})$ sea una matriz definida positiva. La función $q(x)$ es entonces estrictamente convexa y tiene un mínimo único en $x^{(k+1)}$ dado por
$$
\nabla_q(x^{(k+1)})=0
$$
Esto da lugar sistema lineal de ecuaciones:
$$
\nabla f(x^{(k)})=-\nabla^2f(x^{(k)})(x^{(k+1)}-x^{(k)})
$$
La solución para $x_{k+1}$ es el mínimo, lo cual da lugar a la
recurrencia

\begin{equation} \label {eq:newton}
x^{(k+1)}=x^{(k)}-\left[\nabla^2f(x^{(k)})\right]^{-1}\nabla f(x^{(k)})
\end{equation}

En el algoritmo \ref{alg:Newton} se presenta el proceso para encontrar en mínimo de una función utilizando el algoritmo de Newton y abajo la implementación en Matlab

\begin{verbatim}
function r = Newton(g, H, x1)

while 1
    d= inv(H(x1))* g(x1);
    x2 = x1 - d;
    if norm(d) < 0.0001 break;
    else x1 = x2;
    end;
end;
r = x2;
\end{verbatim}

\begin{algorithm}
  \caption{Algoritmo de Newton}
  \label{alg:Newton}
  \begin{algorithmic}
    \REQUIRE{Dado  $f(x)$ y  $x^{(0)}$}
    \ENSURE{$x^{(k)}$}
    \WHILE {$(\| \nabla f(x^{(k)}) \| \leq \varepsilon)$}
     \STATE {$x^{(k+1)} = x^{(k)} - [\nabla^2f(x^{(k)})]^{-1}*\nabla f(x^{(k)})$}
     \STATE {$k \leftarrow k+1$}
    \ENDWHILE \\
    \RETURN{$x^{(k)}$ }
  \end{algorithmic}
\end{algorithm}

\subsection{Ejemplo 1}

Minimizar es $f(x)=x_1e^{-x_1^2-x_2^2}$, con un punto inicial $x_1^{(0)}=-0.5$, $x_2^{(0)}=-0.1$ aplicando el Método de Newton.

El vector gradiente para esta función es:

$$
\nabla f(x_1,x_2) =\left[%
\begin{array}{c}
  (1-2x_1^2)e^{(-x_1^2-x_2^2)} \\
  -2x_1x_2e^{(-x_1^2-x_2^2)} \\
\end{array}
\right]
$$ 

\begin{verbatim}
function g = Gradiente_f1(x) 
n = length(x);
g = zeros(n, 1);

val_e = exp(-x(1)^2-x(2)^2);

g(1) = val_e*(1-2*x(1)^2);
g(2) = val_e*(-2*x(1)*x(2));
\end{verbatim}


El Hessiano es:

$$
\nabla^2 f(x_1,x_2) =
\left[%
\begin{array}{cc}
  (-6x_1+4x_1^3)e^{(-x_1^2-x_2^2)} & (-2x_2+4x_1^2x_2)e^{(-x_1^2-x_2^2)} \\
  (-2x_2+4x_1^2x_2)e^{(-x_1^2-x_2^2)} & (-2x_1 + 4x_1x_2^2)e^{(-x_1^2-x_2^2)} \\
\end{array}
\right]
$$

\begin{verbatim}
function H = Hessiano_f1(x)
n = length(x);
H = zeros(n, n);

val_e = exp(-x(1)^2-x(2)^2);
H(1,1) = val_e*(-6*x(1)+4*x(1)^3);
H(1,2) = val_e*(-2*x(2)+4*x(1)^2*x(2));
H(2,1) = H(1,2);
H(2,2) = val_e*(-2*x(1)+4*x(1)*x(2)^2);
\end{verbatim}

\textbf{Primer iteración}

\begin{equation*}
\left[
\begin{array}{r}
x_1^{(1)}    \\
x_2^{(1)}     
\end{array}
\right] = 
\left[
\begin{array}{r}
   -0.5000 \\
   -0.1000 \\  
\end{array}
\right] -
\left[
\begin{array}{rr}
    1.9276  &  0.0771 \\
    0.0771  &  0.7556 \\
\end{array}
\right]^{-1}
\left[
\begin{array}{r}
    0.3855 \\
   -0.0771
\end{array}
\right] =
\left[
\begin{array}{r}
   -0.7049 \\
    0.0230
\end{array}
\right]
\end{equation*}

\textbf{Segunda iteración}

\begin{equation*}
\left[
\begin{array}{r}
x_1^{(2)}    \\
x_2^{(2)}     
\end{array}
\right] = 
\left[
\begin{array}{r}
   -0.7049 \\
    0.0230
\end{array}
\right] -
\left[
\begin{array}{rr}
    1.7199 &  -0.0002 \\
   -0.0002 &   0.8564\\
\end{array}
\right]^{-1}
\left[
\begin{array}{r}
    0.0038 \\
    0.0197
\end{array}
\right] =
\left[
\begin{array}{r}
   -0.7071 \\
   -0.0000
\end{array}
\right]
\end{equation*}

\textbf{Tercer iteración}

\begin{equation*}
\left[
\begin{array}{r}
x_1^{(3)}    \\
x_2^{(3)}     
\end{array}
\right] = 
\left[
\begin{array}{r}
   -0.7071 \\
   -0.0000
\end{array}
\right] -
\left[
\begin{array}{rr}
    1.7155  &  0.0000\\
    0.0000  &  0.8578
\end{array}
\right]^{-1}
\left[
\begin{array}{r}
    0.0180 \\
   -0.2114
\end{array}
\right] \times 10^{-4} =
\left[
\begin{array}{r}
   -0.7071 \\
   -0.0000
\end{array}
\right]
\end{equation*}

En la siguiente tabla se muestra como la función $f(x)$ disminuye su valor en el proceso de convergencia

\begin{center}
\begin{tabular}{|c|c|}
  \hline
  $k$  & $f(x^{(k)})$ \\
  \hline
   0  & -0.3855  \\
   1  & -0.4287 \\
   2  & -0.4289  \\
   3  & -0.4289  \\
  \hline
\end{tabular}
\end{center}

Para ejecutar dar 

\begin{verbatim}
 
>> Newton(@Gradiente_f1, @Hessiano_f1, [-0.5;-0.1])

ans =

   -0.7071
    0.0000
\end{verbatim}

En la figura \ref{fig:exp} se muestra la función a la derecha y a la izquierda las curvas de nivel. El código en matlab utilizado para generar estas imágenes es

\begin{verbatim}
[x,y] = meshgrid(-4:0.2:4, -4:0.2:4);
z = x.*exp(-x.^2-y.^2);
surf(x,y,z)
\end{verbatim}

\begin{verbatim}
[x,y] = meshgrid(-4:0.2:4, -4:0.2:4);
z = x.*exp(-x.^2-y.^2);
contour(x,y,z)
grid on
\end{verbatim}

\begin{figure}[!htb]
\begin{center}
	\subfigure[Funsión] {\psfig{figure=Capitulo_5/figuras/exp.eps,width=7cm}}
	\subfigure[Curvas de Nivel ] {\psfig{figure=Capitulo_5/figuras/exp_cont.eps,width=7cm}}	
\end{center}
\caption{Mínimo para la función $f(x_1, x_2) = x_1 e^{(-x_1^2 - x_2^2)} $} \label{fig:exp}
\end{figure}
\subsection{Ejemplo 2}

Minimizar es $f(x)=-3.5 x_1 - 2x_2 - x_2^2+x_1^4+2x_1x_2 +x_2^4$, con un punto inicial $x_1^{(0)}=1.0$, $x_2^{(0)}=1.0$ aplicando el Método de Newton.

El vector gradiente para esta función es:

$$
\nabla f(x_1,x_2) =
\left[
\begin{array}{c}
	-3.5 + 4x_1^3 + 2x_2\\
	-2 + 2x_1 - 2x_2 + 4x_2^3\\
\end{array}
\right]
$$ 

\begin{verbatim}
function g = Gradiente_f2(x) 
n = length(x);
g = zeros(n, 1);

g(1) = -3.5 + 4*x(1)^3 + 2*x(2);
g(2) = -2.0 + 2*x(1) - 2*x(2) + 4*x(2)^3;
\end{verbatim}


El Hessiano es:

$$
\nabla^2 f(x_1,x_2) =
\left[
\begin{array}{cc}
  12 x_1^2 & 2\\
  2 & -2+12x_2^2 \\
\end{array}
\right]
$$

\begin{verbatim}
function H = Hessiano_f2(x)
n = length(x);
H = zeros(n, n);

H(1,1) = 12*x(1)^2;
H(1,2) = 2;
H(2,1) = 2;
H(2,2) = -2+12*x(2)^2;
\end{verbatim}

\textbf{Primer iteración}

\begin{equation*}
\left[
\begin{array}{r}
x_1^{(1)}    \\
x_2^{(1)}     
\end{array}
\right] = 
\left[
\begin{array}{r}
     1\\
     1
\end{array}
\right] -
\left[
\begin{array}{rr}
    12  &   2 \\
     2   & 10
\end{array}
\right]^{-1}
\left[
\begin{array}{r}
    2.5000\\
    2.0000
\end{array}
\right] =
\left[
\begin{array}{r}
    0.8190\\
    0.8362
\end{array}
\right]
\end{equation*}

\textbf{Segunda iteración}

\begin{equation*}
\left[
\begin{array}{r}
x_1^{(2)}    \\
x_2^{(2)}     
\end{array}
\right] = 
\left[
\begin{array}{r}
    0.8190 \\
    0.8362
\end{array}
\right] -
\left[
\begin{array}{rr}
    8.0485  &  2.0000\\
    2.0000  &  6.3909
\end{array}
\right]^{-1}
\left[
\begin{array}{r}
    0.3695\\
    0.3044
\end{array}
\right] =
\left[
\begin{array}{r}
    0.7820\\
    0.8001
\end{array}
\right]
\end{equation*}

\textbf{Tercer iteración}

\begin{equation*}
\left[
\begin{array}{r}
x_1^{(3)}    \\
x_2^{(3)}     
\end{array}
\right] = 
\left[
\begin{array}{r}
    0.7820\\
    0.8001
\end{array}
\right] -
\left[
\begin{array}{rr}
    7.3385  &  2.0000\\
    2.0000  &  5.6828
\end{array}
\right]^{-1}
\left[
\begin{array}{r}
    0.0132\\
    0.0129
\end{array}
\right]  =
\left[
\begin{array}{r}
    0.7807\\
    0.7983
\end{array}
\right]
\end{equation*}

\textbf{cuarta iteración}

\begin{equation*}
\left[
\begin{array}{r}
x_1^{(4)}    \\
x_2^{(4)}     
\end{array}
\right] = 
\left[
\begin{array}{r}
    0.7807\\
    0.7983
\end{array}
\right] -
\left[
\begin{array}{rr}
    7.3139  &   2.0000\\
    2.0000  &  5.6483
\end{array}
\right]^{-1}
\left[
\begin{array}{r}
    0.1610\\
    0.3115
\end{array}
\right] \times 10^{-4}  =
\left[
\begin{array}{r}
     0.7807\\
    0.7983
\end{array}
\right] 
\end{equation*}

En la siguiente tabla se muestra como la función $f(x)$ disminuye su valor en el proceso de convergencia

\begin{center}
\begin{tabular}{|c|c|}
  \hline
  $k$  & $f(x^{(k)})$ \\
  \hline
   0  & -2.5000 \\
   1  & -2.9296 \\
   2  & -2.9422  \\
   3  & -2.9423  \\
  \hline
\end{tabular}
\end{center}

Para ejecutar dar 

\begin{verbatim}
 
>> Newton(@Gradiente_f2, @Hessiano_f2, [1;1])

ans =

   -0.7807
    0.7983
\end{verbatim}

En la figura \ref{fig:curvas_eje_02} se muestras las curvas de nivel correspondientes a este ejemplo las cuales fueron determinadas utilizando el siguiente código.

\begin{figure}[!htp]
   \centering
   \centerline{\psfig{figure=Capitulo_5/figuras/fig_eje2.eps,width=10.0cm}}
   \caption{Curvas de nivel para el ejemplo2.}
   \label{fig:curvas_eje_02}
\end{figure}


\begin{verbatim}
[x1, x2] = meshgrid(-2:0.2:2, -2:0.2:2);
f =-3.5*x1-2*x2-x2.^2 + x1.^4 + 2*x1.*x2 + x2.^4;
contour(x1,x2,f)
grid on
xlabel('x_1')
ylabel('x_2')
\end{verbatim}

\subsection{Propiedades del Método de Newton}

Una propiedad importante de este método es que converge en un solo paso cuando se tiene una función estrictamente cuadrática.

\textbf{Demostración:}

Dada una función cuadrática
$$
f(x)=\frac{1}{2}x^TAx + x^Tb + c
$$
el gradiente lo podemos calcular como $\nabla f(x)=Ax+b$ y el Hessiano es $\nabla^2f(x)=A$. Si sustituimos en la ecuación \ref{eq:newton} tenemos
$$
=x^{(k)}-A^{-1}(Ax^{(k)}+b)
$$
$$
x^{(k+1)}=x^{(k)}-x^{(k)}-A^{-1}b
$$
Lo que finalmente da como resultado
$$
x^{(k+1)}=-A^{-1}b
$$

\subsection{Ejemplo 3} 

Dada la función $f(x)= 10x_1^2 +x_2^2$ y un punto inicial $x^{(0)} = [1,2]^T$, calcular el mínimo utilizando el algoritmo de Newton.

Para esta función el gradiente y el Hessiano es:
$$
\nabla f(x) =\begin{array}{|c|}
  \hline
  20x_1 \\
  \hline
  2x_2 \\
  \hline
\end{array}
$$
y el Hessiano es:
$$
A(x)=
\begin{array}{|c|c|}
  \hline
  20 & 0 \\
  \hline
  0 & 2 \\
  \hline
\end{array}%
$$
Aplicando la ecuación \ref{eq:newton} tenemos
$$
x^{(1)}=
\begin{array}{|c|}
\hline
1 \\
\hline
2 \\
\hline
\end{array}
-
\begin{array}{|c|c|}
  \hline
  20 & 0 \\
  \hline
  0 & 2 \\
  \hline
\end{array}%
^{-1}
\begin{array}{|c|}
  \hline
  20 \\
  \hline
  4 \\
  \hline
\end{array} =
\begin{array}{|c|}
  \hline
  0 \\
  \hline
  0 \\
  \hline
\end{array}
$$

\subsection{Ejemplo 4}

Minimizar es $f(x)=x_1e^{-x_1^2-x_2^2}$, con un punto inicial $x_1^{(0)}=-0.5$, $x_2^{(0)}=-0.5$ aplicando el Método de Newton.

\textbf{Primer iteración}

\begin{equation*}
\left[
\begin{array}{r}
x_1^{(1)}    \\
x_2^{(1)}     
\end{array}
\right] = 
\left[
\begin{array}{r}
   -0.5000 \\
   -0.5000
\end{array}
\right] -
\left[
\begin{array}{rr}
    1.5163  &   0.3033 \\
    0.3033  &  0.3033
\end{array}
\right]^{-1}
\left[
\begin{array}{r}
    0.3033 \\
   -0.3033
\end{array}
\right] =
\left[
\begin{array}{r}
   -1.0000 \\
    1.0000
\end{array}
\right]
\end{equation*}

\textbf{Segunda iteración}

\begin{equation*}
\left[
\begin{array}{r}
x_1^{(2)}    \\
x_2^{(2)}     
\end{array}
\right] = 
\left[
\begin{array}{r}
   -1.0000 \\
    1.0000
\end{array}
\right] -
\left[
\begin{array}{rr}
    0.2707  &  0.2707\\
    0.2707  & -0.2707
\end{array}
\right]^{-1}
\left[
\begin{array}{r}
   -0.1353 \\
    0.2707
\end{array}
\right] =
\left[
\begin{array}{r}
   -1.2500 \\
    1.7500
\end{array}
\right]
\end{equation*}

\textbf{Tercer iteración}

\begin{equation*}
\left[
\begin{array}{r}
x_1^{(3)}    \\
x_2^{(3)}     
\end{array}
\right] = 
\left[
\begin{array}{r}
   -1.2500 \\
    1.7500
\end{array}
\right] -
\left[
\begin{array}{rr}
   -0.0031  &  0.0729 \\
    0.0729  & -0.1256
\end{array}
\right]^{-1}
\left[
\begin{array}{r}
   -0.0208 \\
    0.0429
\end{array}
\right]  =
\left[
\begin{array}{r}
   -1.3535 \\
    2.0314
\end{array}
\right]
\end{equation*}

\textbf{Cuarta iteración}

\begin{equation*}
\left[
\begin{array}{r}
x_1^{(4)}    \\
x_2^{(4)}     
\end{array}
\right] = 
\left[
\begin{array}{r}
   -1.3535\\
    2.0314
\end{array}
\right] -
\left[
\begin{array}{rr}
   -0.0046  &  0.0280\\
    0.0280  & -0.0507
\end{array}
\right]^{-1}
\left[
\begin{array}{r}
   -0.0069\\
    0.0142
\end{array}
\right]  =
\left[
\begin{array}{r}
   -1.4416\\
    2.2629
\end{array}
\right]
\end{equation*}


Para ejecutar dar 

\begin{verbatim}
 
>> Newton(@Gradiente_f1, @Hessiano_f1, [-0.5;-0.5])

\end{verbatim}


En la siguiente tabla se muestra como la función $f(x)$ disminuye su valor en el proceso de convergencia

\begin{center}
\begin{tabular}{|r|r|r|r|}
  \hline
  % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
  $k$  & $x_1^{(k)}$ & $x_2^{(k)}$ & $f(x^{(k)})$ \\
  \hline
    0  &  -0.5000  &  -0.5000 &  -0.3033\\
    1  & -1.0000   & 1.0000   & -0.1353\\
    2  & -1.2500   & 1.7500   & -0.0123\\
    3  & -1.3535   &  2.0314  &  -0.0035\\
    4  & -1.4416   &  2.2629  &  -0.0011\\
    5  & -1.5206   &  2.4654  &  -0.0003\\
    6  & -1.5933   &  2.6481  & -0.0001\\
    7  & -1.6611   &  2.8161  & -0.0000\\
    8  & -1.7251   &  2.9727  &  -0.0000\\
    9  & -1.7859   &  3.1198  &  -0.0000\\
  \hline
\end{tabular}
\end{center}
De los datos de la tabla podemos notar que la función en lugar de disminuir crece. Esta condición se debe a que el punto inicial para el método de Newton esta muy lejos y la aproximación cuadrática, resulta in apropiada. En la figura \ref{fig:xexp} se muestra la función y se puede ver que esta tiene solamente un mínimo.

\begin{figure}[!htp]
   \centering
   \centerline{\psfig{figure=Capitulo_5/figuras/fun_xexp.eps,width=8.0cm}}
   \caption{Función $f(x)=x_1e^{-x_1^2-x_2^2}$.}
   \label{fig:xexp}
\end{figure}


\subsection{Problemas de convergencia del Método de Newton}

Consideremos el caso de una función donde se tienen múltiples raíces como es el caso de la función $f(x) = (x-3)^5/5$ cuyo primer derivada es $f'(x) = (x-3)^4$. Note que en este caso la solución es un polinomio de grado 4 en cuyo caso tendrá cuatro soluciones con $x^* = 3$. Si aplicamos el método de Newton tendremos una convergencia como se muestra a continuación
 
\begin{center}
\begin{tabular}{r|r}
    $k$ & $f(x^{(k)})$ \\
    \hline
    1 &    2.2500 \\
    2 &    2.4375 \\
    3 &    2.5781 \\
    4 &    2.6836 \\
    5 &    2.7627 \\
    6 &    2.8220 \\
    7 &    2.8665 \\
    8 &    2.8999 \\
    9 &    2.9249 \\
   10 &    2.9437 \\
   \vdots & \vdots \\
    44 &    3.0000\\
    45 &    3.0000\\
\end{tabular}
\end{center}

Si uno desea aplicar el método a una función arbitraria encontraremos algunas dificultades, esencialmente debido al hecho de que no encontraremos convergencia global. Si el punto inicial $x^{(0)}$ esta muy lejos de $x^*$ el método no podrá converger.


Para comenzar, la aproximación de $f(x$) dada por $q(x)$ es solamente valida a la vecindad de $x_k$, el tamaño de paso puede ser controlado a través de una formula iterativa de tipo:

\begin{center}
$x^{(k+1)}=x^{(k)}-\lambda^{(k)}\left[\nabla^2f(x^{(k)})\right]^{-1}\nabla
f(x^{(k)})$
\end{center}

donde $\lambda^{(k)}$ es un escalar, seleccionado por ejemplo, de tal manera que $||x^{(k+1)}-x^{(k)}||$ no sea muy largo. También podemos seleccionarlo tal que $\lambda^{(k)}$ minimice $\phi(\lambda^{(k)})=f(x^{(k)}+\lambda^{(k)}d^{(k)})$ en la dirección de:

\begin{center}
$d^{(k)}=-\left[\nabla^2f(x^{(k)})\right]^{-1}\nabla f(x^{(k)})$
\end{center}

\subsection{Ejemplo 5}

Minimizar $f(x)=\frac{(x-3)^5}{5}$ utilizando el método de Newton y la modificación planteada utilizando un valor $\lambda^{(k)} = 2.9$. Considere un punto inicial $x^{(0)} = 2$.
\begin{center}
\begin{tabular}{|c|c|c|}
  \hline
  % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
  k & Newton & Newton con $\lambda=2.9$ \\
  \hline
  0 &  2.00000 & 2.00000 \\
  1 &  2.33333 & 2.96667 \\
  2 &  2.55556 & 2.99889 \\
  3 &  2.70370 & 2.99996 \\
  4 &  2.80247 & 3.00000 \\
  5 &  2.86831 & 3.00000 \\
  6 &  2.91221 & 3.00000 \\
  7 &  2.94147 & 3.00000 \\
  8 &  2.96098 & 3.00000 \\
  9 &  2.97399 & 3.00000 \\
 10 &  2.98266 & 3.00000 \\
 11 &  2.98844 & 3.00000\\
  \hline
\end{tabular}
\end{center}

\subsection{Método de Newton Modificado}

Una alternativa propuesta por Ralston y Rabinowitz (1978) es la de definir una nueva función $u(x)$, que es el cociente del gradiente de $g(x)$ y su derivada (en una dimensión)
$$
u'(x) = \frac{f'(x)}{f''(x)}
$$
Se puede mostrar que esta función tiene las mismas raíces que $f'(x)$ y que la multiplicidad de raíces no afectará. La
formulación del método de Newton es:
$$
x_{k+1} = x_k - \frac{u'(x)}{u''(x)}
$$
La derivada de $u'(x)$ es:
$$
u''(x) = \frac{f''(x)f''(x) - f'(x)f'''(x)}{[f''(x)]^2}
$$
Sustituyendo esta, tenemos la formulación final del Método de
Newton modificado.
$$
x^{(k+1)} = x^{(k)} - \frac{f''(x^{(k)})f'(x^{(k)})}{f''(x^{(k)})f''(x^{(k)}) - f'(x^{(k)})f'''(x^{(k)})}
$$

\begin{verbatim}
function y = Newton_1d_modificado(der1f, der2f, der3f, x)
x1 = x;
x2 = x;
for iter = 1:10
    delta = der1f(x1)*der2f(x1)/(der2f(x1)^2 - der1f(x1)*der3f(x1));
    x2 = x1 - delta;
    disp([iter x2]);
    if abs(delta) < 1e-6 break;
    else
        x1 =x2;
    end;
end;
y = x2;
\end{verbatim}

\subsubsection{Ejemplo}

Para la función $f(x) = \frac{(x-3)^5}{5}$ para un valor inicial $x_0 =2$ la solución es:

La primer derivada es:
$$
f'(x) = (x-3)^4
$$
La segunda derivada
$$
f''(x) = 4(x-3)^3
$$
y la tercer derivada 
$$
f'''(x)= 12(x-3)^2
$$

Evaluando en el punto inicial tenemos $f'(2) = 1$, $f''(2) = -4$ y $f'''(2) = 12$ sustituyendo en la iteración de Newton tenemos

$$
x^{(1)} = 2 - \frac{(1)\times (-4)}{(4)^2 - (1)*(12)} = 3
$$
Note que la solución se da en una sola iteración. Para ejecutar hacer


\begin{verbatim}
Newton_1d_modificado(@der_fun_3, @der2_fun_3, @der3_fun_3,2)
\end{verbatim}

\section[Metodo Simplex]{Optimización lineal restringida (programación lineal) Método Simplex}

En programación lineal (o Pl por simplicidad) es un procedimiento de optimización que trata de cumplir con un objetivo como maximizar las utilizares o minimizar el costo, en presencia de restricciones como lo son las fuentes limitadas. El término lineal denota que las funciones matemáticas que representan ambos, la función objetivo y las restricciones son lineales. 

\subsection{Forma estándar}

El problema básico de programación lineal consiste en dos partes principales: la función objetivo y un conjunto de restricciones. Para un problema de maximización la función objetivo es por lo general expresada como:

Maximizar

$$
Z = c_1 x_1 + c_2 x_2 + \cdots + c_n x_n
$$

donde $c_j$ es el pago por cada unidad de la $j$-esima actividad que se lleva a cabo y $x_j$ es la magnitud o cantidad de la $j$-esima actividad.

Las restricciones se pueden representar en forma general como:

$$
a_{i,1} x_1 + a_{i,2} x_2 + \cdots + a_{i,n} x_n\leq b_i
$$

donde $a_{i,j}$ cantidad de $i$-ésima fuente que se consume por cada unidad de $j$-ésima actividad y $b_i$ cantidad de la $i$-ésima fuente que ésta disponible.

Suponga que una planta procesadora de gasolina recibe cada semana una cantidad fija de materia prima. Esta última se procesa en dos tipos de gasolina de calidad regular y premium. Estas clases de gasolina son de alta demanda; es decir, se tiene garantizada su venta y se obtiene diferentes utilidades para la compañía. Si embargo su producción involucra ambas restricciones, tiempo y almacenaje en sitio. Por ejemplo, sólo una de las clases se puede producir a la vez, y las instalaciones están abiertas solo 80 horas. Además, existe un limite de almacenamiento para cada uno de los productos. Todos estos factores se en listan en la siguiente tabla:


\begin{center}
\begin{tabular}{|c|c|c|c|}
\hline
Recurso & Regular & Premium & Disponibilidad \\
\hline
\hline
Materia Prima & 7 m$^3$/tonelada & 11 m$^3$/tonelada & 77 m$^3$/semana\\
Tiempo de Producción & 10 hr/tonelada & 8 hr/tonelada & 80 hr/semana\\
Almacenamiento & 9 toneladas & 6 toneladas & \\
\hline
\hline
Aprovechamiento & 150 pesos/tonelada & 175 pesos/tonelada & \\
\hline
\end{tabular}
\end{center}
 
 \textbf{Solución}
 
El ingeniero  que opera esta planta debe decidir la cantidad a producir de cada gasolina para maximizar las utilidades. Si las cantidades producidas cada semana de gasolina regular y premium son designadas $x_1$ y $x_2$ respectivamente, la ganancia total se puede calcular como
 
Ganancia total  $=150 x_1 + 175 x_2$
 
o escribirla como una función objetivo en programación lineal
 
Maximiza $Z = 150 x_1 + 175 x_2$
 
Las restricciones se pueden desarrollar en una forma similar. Por ejemplo el total de gasolina cruda utilizada se puede calcular como
 
Total de materia prima  utilizada $ =  7x_1 + 11 x_2$ 

Este total no puede exceder el abastecimiento disponible de 77 m3/semana, así que la restricción se puede representar como

$$
7x_1 + 11 x_2 \leq 77
$$

Las restricciones restantes se pueden desarrollar en una forma similar, la formulación total resultante para la PL está dada por

Maximizar $Z = 150 x_1 + 175 x_2$

Sujeta a 

\begin{eqnarray*}
7 x_1 + 11 x_2 \leq 77\\
10 x_1 + 8 x_2 \leq 80\\
x_1 \leq 9\\
x_2 \leq 6\\
x_1 \geq 0\\
x_2 \geq 0
\end{eqnarray*}

\subsection{Método Simplex}

El método Simplex se basa en la suposición de que la solución óptima estará en un punto extremo. Así, el procedimiento debe debe de ser capaz de discernir si durante la solución a un problema ocurre un punto extremo. Para realizar esto, las ecuaciones con restricciones se reformulan como igualdades por medio de la introducción de las llamadas variables de holgura. 

Variables de Holgura. Como lo indica su nombre, una variable de holgura mide cuanto de una fuente restringida esta disponible es decir, cuanto holgura esta disponible. Por ejemplo, recuerde la fuente restringida que se utilizo en el ejemplo: 

$$
7 x_1 + 11 x_2 \leq 77
$$

Se puede definir una variable $h_1$ como la cantidad de gasolina cruda que no se usa para una nivel de producción en particular $(x_1, x_2)$. Si esta cantidad se agrega al lado izquierdo de la restricción, forma la relación exacta. 

$$
7 x_1 + 11 x_2 + h_1 = 77
$$

Ahora se reconoce lo que la variable de holgura nos indicaba. Si esta es positiva, significa que tiene algo de holgura para esta restricción. Esto es, se cuenta con algo más de recursos que no han sido utilizado por completo. Si es negativa nos indica que nos hemos excedido en la restricción. Finalmente, si es cero denota que se cumplió exactamente con la restricción. Es decir, se dispuso de todo el recurso puesto que esta es la condición donde las líneas de restricción se interceptan, la variable de holgura proporciona un medio para determinar puntos extremos.

Una variable de holgura diferente se desarrolla para cada ecuación restringida resultando en lo que se llama versión completamente aumentada.

$$
Z = 150 x_1 + 175 x_2
$$

Sujeta a 
\[
\begin{array}{rrrrrrr}
    7 x_1 & + 11 x_2  & + h_1 &           &           &          & = 77   \\
  10 x_1 &  + 8 x_2  &            & + h_2 &           &          &=80 \\
         x_1&               &            &           & + h_3 &          &= 9\\
               &  x_2       &            &           &            & + h_4 & = 6
\end{array}
\]

Advierta como se han formulado de igual modo las cuatro ecuaciones, de tal manera que las incógnitas están alineadas en las columnas. Se hizo así para resaltar que ahora se trata de un sistema de ecuaciones algebraicas lineales.


Los pasos del Método Simplex son los siguientes:

\begin{enumerate}
\item Utilizando la forma estándar, determinar una solución básica factible inicial igualando a las $n-m$ variables igual a cero (el origen).
\item  Seleccionar la variable de entrada de las variables no básicas que al incrementar su valor pueda mejorar el valor en la función objetivo. Cuando no exista esta situación, la solución actual es la óptima; si no, ir al siguiente paso. Esta selección se hada tomando el $j$-ésimo coeficiente de la función objetivo con el coeficiente $c_j$ mas grande.
\item Seleccionar la variable de salida de las variables básicas actuales. Tomamos la $i$-ésima variable tal que $b_i/a_{i,j}$ tenga el coeficiente mas pequeño no negativo.
\item  Determinar la nueva solución al hacer la variable de entrada básica y la variable de salida no básica, ir al paso 2 (actualizar). La solución se obtiene haciendo una iteración de Gauss-Jordan.
\end{enumerate}

La implementación en Matlab es:

\begin{verbatim}
function x = Simplex(C, A, b)
N = length(b);
B  = [A,  diag(ones(N,1)), b];
B = [B; C, zeros(1, N+1)]
M = length(B(1,:));

while 1 
    [mm, m] = max(B(N+1,:)) ;
    [mn, n] = imin(B(1:N, M)./B(1:N, m));
    
    if mm == 0
        break;
    end;
      
    B(n,:) = B(n, :)./B(n,m);
    
    for k = 1: N+1
        if n ~= k
            B(k,:) = B(k,:) - B(n,:)*B(k,m);
        end;
    end; 
end;
x=B;\end{verbatim}

\subsection{Ejemplo 1}
Mostrar la región de factibilidad del  problema de la Gasolina planteado y resolver utilizando el método de Simplex.  

$$
Z = 150 x_1 + 175 x_2
$$

Sujeta a 
\[
\begin{array}{rrrrrrr}
    7 x_1 & + 11 x_2  & + h_1 &           &           &          & = 77   \\
  10 x_1 &  + 8 x_2  &            & + h_2 &           &          &=80 \\
         x_1&               &            &           & + h_3 &          &= 9\\
               &  x_2       &            &           &            & + h_4 & = 6
\end{array}
\]

Para determinar la región de factibilidad se corrió el siguiente código en matlab y en la Figura \ref{fig:eje_01} se muestra en verde la región de factibilidad, las restricciones en negro y en rojo la función objetivo con $Z = 1000$.

\begin{verbatim}
	hold off;

	[x1,x2] = meshgrid(-1:0.1:10, -1:0.1:10);
	f1 = 7*x1 + 11*x2 - 77;
	f2 = 10*x1 + 8*x2 - 80;
	f3 = x1 - 9;
	f4 = x2 - 6;
	f5 = x1;
	f6 = x2;

	Z = 150*x1 + 175*x2 - 1000;

	contour(x1, x2, f1, [0,0], 'k');
	hold on;
	contour(x1, x2, f2, [0,0], 'k');
	contour(x1, x2, f3, [0,0], 'k');
	contour(x1, x2, f4, [0,0], 'k');
	contour(x1, x2, f5, [0,0], 'k');
	contour(x1, x2, f6, [0,0], 'k');
	contour(x1, x2, Z, [0,0], 'r');
	xlabel('x_1');
	ylabel('x_2');
\end{verbatim}

\begin{figure}[!htp]
   \centering
   \centerline{\psfig{figure=Capitulo_5/figuras/eje_01.jpg,width=14.0cm}}
   \caption{Para el ejemplo 1, en verde área de Factibilidad con la función objetivo en rojo para Z = 100. }
   \label{fig:eje_01}
\end{figure}

Comenzamos por escribir la forma estándar del problema seleccionando como solución básica factible $x_1 = 0$, $x_2 =0$, $h_1 =77$, $h_2 =80$, $h_3 =9$ y $h_4 =6$. Note que las variables de holgura tiene los valores que hacen las desigualdades ver como igualdades

\[
\left[
\begin{array}{rr|rrrr|r}
           7.00      &    11.00     &      1     &       0      &       0      &       0    &     77.00\\
         10.00      &      8.00     &      0     &       1      &       0      &       0    &     80.00\\
           1.00      &       0         &      0     &        0     &       1      &       0    &       9.00\\
                0      &      1.00     &      0     &        0     &       0      &       1    &       6.00\\
         \hline
       150.00      &  175.00     &      0     &        0     &       0      &       0    &       0\\
\end{array}
\right]
\]

Note que en nuestra formulación estándar la función objetivo aparece con los coeficientes negativos, por esta razón repetiremos el procedimiento hasta que no tengamos ningún término negativo en el renglón correspondiente a la función de costo.

\textbf{Primer iteración}

Dado que el costo mayor es 175 comenzamos por introducir la variable $x_2$  por lo tanto la columna es la 2 y sacar la variable de holgura $h_4$ dado que esta restricción es la mas cercana del punto básico. Por lo tanto el pivote es 4, 2
          
\[
\left[
\begin{array}{r>{\columncolor{gray}}r|rrrr|r}
           7.00      &     11.00      &      1     &       0      &       0      &       0    &     77.00\\
         10.00      &       8.00      &      0     &       1      &       0      &       0    &     80.00\\
           1.00      &        0          &      0     &        0     &       1      &       0    &       9.00\\
           \rowcolor{gray} 
                0      &       1.00      &      0     &        0     &       0      &      1      &       6.00\\
         \hline
       150.00      &   175.00      &      0     &        0     &       0      &       0    &       0\\
\end{array}
\right]
\begin{array} {r}
7\\
10\\
\infty \\
6 \\
\\
\end{array}
\]
 
Una vez aplicado el método de Gauss Jordan tenemos
          
\[
\left[
\begin{array}{rr|rrrr|r}
          7.00       &      0     &     1.00    &         0     &        0    &    -11.00   &      11.00\\
         10.00      &       0     &        0     &     1.00    &         0   &      -8.00   &      32.00\\
          1.00       &      0      &       0      &       0       &   1.00    &         0       &   9.00\\
          0         &   1.00     &      0      &      0       &      0      &   1.00     &    6.00\\
             \hline
        150.00     &        0    &         0    &         0     &        0    &    -175.00  &     -1050.00\\
\end{array}
\right]
\]

\textbf{Segunda Iteración}

La siguiente variable a introducir es $x_1$ dado que el mayor de los costos es 150 y sacar la variable de holgura $h_1$  dado que es la mas cercana de nuestro punto extremo. Por lo tanto tenemos el elemento 1, 1 como pivote

\[
\left[
\begin{array}{>{\columncolor{gray}}rr|rrrr|r}
	    \rowcolor{gray} 
          7.00       &      0     &     1.00    &         0     &        0    &    -11.00   &      11.00\\
         10.00      &       0     &        0     &     1.00    &         0   &      -8.00   &      32.00\\
          1.00       &      0      &       0      &       0       &   1.00    &         0       &   9.00\\
          0         &   1.00     &      0      &      0       &      0      &   1.00     &    6.00\\
             \hline
        150.00     &        0    &         0    &         0     &        0    &    -175.00  &     -1050.00\\
\end{array}
\right]
\begin{array} {r}
1.57\\
3.2\\
9 \\
\infty \\
\\
\end{array}
\]

Aplicando el método de Gauss-Jordan en el elemento 1,1 tenemos
\[
\left[
\begin{array}{rr|rrrr|r}
          1.00     &        0    &      0.14    &         0     &        0    &     -1.57      &    1.57\\
             0       &      0      &   -1.43      &    1.00     &         0    &      7.71    &     16.29\\
             0       &      0      &   -0.14      &       0       &   1.00     &     1.57     &     7.43\\
             0       &   1.00    &         0       &      0        &     0       &   1.00       &   6.00\\
             \hline
             0       &      0      &   -21.43      &       0       &      0     &    60.71    &  -1285.71\\
\end{array}
\right]
\]

\textbf{Tercer iteración}

En la función de Costo tenemos un término positivo de 60.71 corresponden a la variable de holgura $h_4$, por lo tanto introducimos la variable $h_4$ y sacamos la variable $h_2$. Por lo tanto el pivote para esta iteración es  2,  6

\[
\left[
\begin{array}{rr|rrr>{\columncolor{gray}}r|r}
          1.00     &        0    &      0.14    &         0     &        0    &     -1.57      &    1.57\\
           \rowcolor{gray} 
             0       &      0      &   -1.43      &    1.00     &         0    &      7.71    &     16.29\\
             0       &      0      &   -0.14      &       0       &   1.00     &     1.57     &     7.43\\
             0       &   1.00    &         0       &      0        &     0       &   1.00       &   6.00\\
             \hline
             0       &      0      &  -21.43      &       0       &      0     &    60.71    &   -1285.71\\
\end{array}
\right]
\begin{array} {r}
-1.00\\
2.11\\
4.73 \\
6 \\
\\
\end{array}
\]

Una vez eliminada la variable de holgura $h_2$ tenemos

\[
\left[
\begin{array}{rr|rrrr|r}
          1.00     &        0     &    -0.15     &     0.20     &        0     &        0     &     4.89\\
             0       &      0       &  -0.19       &   0.13       &      0       &   1.00     &     2.11\\
             0       &      0       &   0.15       &  -0.20       &   1.00     &        0     &     4.11\\
             0       &   1.00     &     0.19     &    -0.13     &        0     &        0     &     3.89\\
             \hline
             0       &      0       & -10.19      &   -7.87       &      0       &      0     & -1413.89
\end{array}
\right]
\]

Dado que no tenemos ningún termino de costo negativo, la solución la podemos obtener de la formulación de la siguiente manera. Para $x_1$ la solución esta en el renglón 1 dado que hay un 1 en la posición 1,1 y es $x_1 = 4.89$. Para $x_2$ el 1 esta en la posición 4 por lo tanto la solución es $x_2 = 3.89$. Para $h_1$ no tenemos un 1 en la columna correspondiente y de igual manera para $h_2$, por lo tanto tenemos $h_1 =0$ y $h_2 =0$. Para $h_3$ tenemos un 1 en el renglón 3 por lo tanto $h_3 = 4.11$. Finalmente para $h_4$ tenemos un 1 en la posición 2 por lo tanto $h_4= 2.11$. El costo total es $Z = 1413.89$ y lo podemos verificar haciendo $Z = 150 \times 4.89 + 175 \times 3.89 = 1413.89$. La solución en Matlab se calcula como

\begin{verbatim}
>> B = Simplex([150 175], [7 11; 10 8; 1 0; 0 1], [77; 80; 9; 6])

B =

      7.00     11.00      1.00         0         0         0     77.00
     10.00      8.00         0      1.00         0         0     80.00
      1.00         0         0         0      1.00         0      9.00
         0      1.00         0         0         0      1.00      6.00
    150.00    175.00         0         0         0         0         0


B =

      1.00         0     -0.15      0.20         0         0      4.89
         0         0     -0.19      0.13         0      1.00      2.11
         0         0      0.15     -0.20      1.00         0      4.11
         0      1.00      0.19     -0.13         0         0      3.89
         0         0    -10.19     -7.87         0         0  -1413.89
\end{verbatim}

\subsection{Ejemplo 2}

Maximizar la función

\begin{eqnarray*}
Z = 6 x_1 + 4 x_2\\
\texttt{Sujeta a }\\
2 x_1 + 2 x_2 \leq 160\\
x_1 + 2x_2 \leq 120\\
4 x_1 + 2 x_2  \leq 280
\end{eqnarray*}

En la Figura \ref{fig:eje_02}, se muestra en verde la región de factibilidad y en rojo la función objetivo para $Z = 350$. El código para trazar la figura se muestra a continuacion:

\begin{verbatim}
hold off;

[x1,x2] = meshgrid(-10:1:120, -10:1:120);
f1 = 2*x1 + 2*x2 - 160;
f2 = x1 + 2*x2 -120;
f3 = 4*x1 + 2*x2 -280;
f4 = x1;
f5 = x2;

Z = 6*x1 + 4*x2 - 350;

contour(x1, x2, f1, [0,0], 'k');
hold on;
contour(x1, x2, f2, [0,0], 'k');
contour(x1, x2, f3, [0,0], 'k');
contour(x1, x2, f4, [0,0], 'k');
contour(x1, x2, f5, [0,0], 'k');
contour(x1, x2, Z, [0,0], 'r');
xlabel('x_1');
ylabel('x_2');

\end{verbatim}

\begin{figure}[!htp]
   \centering
   \centerline{\psfig{figure=Capitulo_5/figuras/eje_02.jpg,width=14.0cm}}
   \caption{Ejemplo 2, Área de Factibilidad en verde y en rojo la función objetivo con Z=350.}
   \label{fig:eje_02}
\end{figure}

Su formulación estándar es:

\[
\left[
\begin{array}{rr|rrr|r}
          2.00     &     2.00     &     1.00     &        0       &      0     &   160.00\\
          1.00     &     2.00     &        0       &   1.00       &      0     &   120.00\\
          4.00     &     2.00     &        0       &      0         & 1.00     &   280.00\\
          \hline
          6.00     &     4.00     &        0       &      0         &    0       &      0\\
\end{array}
\right]
\]

\textbf{Primer iteración}

Tomamos como pivote el renglón 3 y columna 1
\[
\left[
\begin{array}{>{\columncolor{gray}}rr|rrr|r}
          2.00     &     2.00     &     1.00     &        0       &      0     &   160.00\\
          1.00     &     2.00     &        0       &   1.00       &      0     &   120.00\\
          \rowcolor{gray}
          4.00     &     2.00     &        0       &      0         & 1.00     &   280.00\\
          \hline
          6.00     &     4.00     &        0       &      0         &    0       &      0\\
\end{array}
\right]
\begin{array} {r}
80.0\\
120.0\\
70 \\
\\
\end{array}
\]

Realizada la eliminación de la variable $x_1$ tenemos

\[
\left[
\begin{array}{rr|rrr|r}
             0     &     1.00     &     1.00     &        0      &   -0.50    &     20.00\\
             0     &     1.50     &        0       &   1.00      &   -0.25    &     50.00\\
          1.00   &       0.50   &          0     &        0      &    0.25     &    70.00\\
          \hline
             0     &     1.00     &        0       &      0        & -1.50      & -420.00\\
\end{array}
\right]
\]

\textbf{Segunda Iteración}

Seleccionamos como pivote al renglón 1 columna 2.
\[
\left[
\begin{array}{r>{\columncolor{gray}}r|rrr|r}
            \rowcolor{gray}
             0     &     1.00     &     1.00     &        0      &   -0.50    &     20.00\\
             0     &     1.50     &        0       &   1.00      &   -0.25    &     50.00\\
          1.00   &       0.50   &          0     &        0      &    0.25     &    70.00\\
          \hline
             0     &     1.00     &        0       &      0        &  -1.50      &  -420.00\\
\end{array}
\right]
\begin{array} {r}
20.00\\
33.33\\
70.00 \\
\\
\end{array}
\]
Realizada la eliminación de $x_2$ tenemos          

\[
\left[
\begin{array}{rr|rrr|r}
             0    &      1.00     &     1.00      &       0     &    -0.50    &     20.00\\
             0     &        0       &  -1.50        &  1.00     &     0.50    &     20.00\\
          1.00   &          0     &    -0.50      &       0     &     0.50    &     60.00\\
          \hline
             0     &        0       &   -1.00         &    0       &  -1.00      &  -440.00\\
\end{array}
\right]
\]

Dado que no tenemos términos negativos en el renglón de costos terminamos con solución   $x_1 = 60$, $x_2 = 20$, $h_1 =0$, $h_2 =20$, $h_3 =0$ y costo $Z = 440$ 

Para resolver en Matlab hacer

\begin{verbatim}
>> Simplex([6 4], [2 2; 1 2; 4 2], [160; 120; 280])

B =

     2     2     1     0     0   160
     1     2     0     1     0   120
     4     2     0     0     1   280
     6     4     0     0     0     0


ans =

         0    1.0000    1.0000         0   -0.5000   20.0000
         0         0   -1.5000    1.0000    0.5000   20.0000
    1.0000         0   -0.5000         0    0.5000   60.0000
         0         0   -1.0000         0   -1.0000 -440.0000

>> 
\end{verbatim}

\subsection{Ejemplo 3}

Maximizar la función

\begin{eqnarray*}
Z = 3 x_1 + 4 x_2\\
\texttt{Sujeta a }\\
2 x_1 + 3 x_2 \leq 1200\\
2 x_1 + x_2 \leq 1000\\
4 x_2  \leq 800
\end{eqnarray*}

La formulación utilizando Simplex queda

\[
\left[
\begin{array}{rr|rrr|r}
          2.00    &     3.00      &    1.00      &       0       &      0    &   1200.00\\
          2.00    &      1.00     &        0       &   1.00      &       0   &    1000.00\\
             0      &    4.00       &      0         &    0        &  1.00     &   800.00\\
             \hline
          3.00    &      4.00     &        0       &      0      &       0     &        0
\end{array}
\right]
\]

\textbf{Primer Iteración}

Tomamos como pivote el renglón 3 y columna 2
\[
\left[
\begin{array}{r>{\columncolor{gray}}r|rrr|r}
          2.00    &     3.00      &    1.00      &       0       &      0    &   1200.00\\
          2.00    &      1.00     &        0       &   1.00      &       0   &    1000.00\\
          \rowcolor{gray}
             0      &    4.00       &      0         &    0        &  1.00     &   800.00\\
             \hline
          3.00    &      4.00     &        0       &      0      &       0     &        0
\end{array}
\right]
\begin{array} {r}
400.00\\
1000.00\\
200.00 \\
\\
\end{array}
\]

Reduciendo tememos como resultado

\[
\left[
\begin{array}{rr|rrr|r}
          2.00     &        0     &     1.00     &        0     &    -0.75    &   600.00\\
          2.00     &        0     &        0       &   1.00     &    -0.25    &    800.00\\
             0       &   1.00     &        0       &      0      &    0.25      &  200.00\\
             \hline
          3.00     &        0     &        0       &      0      &    -1.00     &   -800.00\\
\end{array}
\right]
\]

\textbf{Segunda Iteración }

Tomamos como pivote el renglón 1 y columna  1

\[
\left[
\begin{array}{>{\columncolor{gray}}rr|rrr|r}
	   \rowcolor{gray}
          2.00     &        0     &     1.00     &        0     &    -0.75    &   600.00\\
          2.00     &        0     &        0       &   1.00     &    -0.25    &    800.00\\
             0       &   1.00     &        0       &      0      &    0.25      &  200.00\\
             \hline
          3.00     &        0     &        0       &      0      &    -1.00     &   -800.00\\
\end{array}
\right]
\begin{array} {r}
300.00\\
400.00\\
\infty \\
\\
\end{array}
\]

Aplicando la eliminación tenemos

\[
\left[
\begin{array}{rr|rrr|r}
          1.00      &       0      &    0.50     &        0     &    -0.38    &    300.00\\
             0        &     0        & -1.00       &   1.00     &     0.50    &    200.00\\
             0        &  1.00      &       0        &     0        &  0.25      &  200.00\\
             \hline
             0        &     0        &  -1.50        &     0       &   0.12      & -1700.00\\
\end{array}
\right]
\]

\textbf{Tercer iteración}

Tomamos como pivote el renglón 2 la columna  5
\[
\left[
\begin{array}{rr|rr>{\columncolor{gray}}r|r}
          1.00      &       0      &    0.50     &        0     &    -0.38    &    300.00\\
             \rowcolor{gray}
             0        &     0        & -1.00       &   1.00     &     0.50    &    200.00\\
             0        &  1.00      &       0        &     0        &  0.25      &  200.00\\
             \hline
             0        &     0        &  -1.50        &     0       &  0.12      & -1700.00\\
\end{array}
\right]
\begin{array} {r}
-789.47\\
400.00\\
800.00 \\
\\
\end{array}
\]

\[
\left[
\begin{array}{rr|rrr|r}
          1.00      &       0   &      -0.25    &      0.75     &        0    &    450.00\\
             0        &     0     &    -2.00      &    2.00       &   1.00    &    400.00\\
             0        &  1.00   &       0.50    &     -0.50     &        0    &    100.00\\
             \hline
             0        &     0     &     -1.25      &    -0.25       &      0      & -1750.00
\end{array}
\right]
\]

La solución es $x_1 = 450$, $x_2 = 100$ con un costo $Z = 1750.00$

La solución en Matlab se calcula haciendo

\begin{verbatim}

>> B = Simplex([3 4], [2 3; 2 1; 0 4], [1200; 1000; 800])

B =

      2.00      3.00      1.00         0         0   1200.00
      2.00      1.00         0      1.00         0   1000.00
         0      4.00         0         0      1.00    800.00
      3.00      4.00         0         0         0         0


B =

      1.00         0     -0.25      0.75         0    450.00
         0         0     -2.00      2.00      1.00    400.00
         0      1.00      0.50     -0.50         0    100.00
         0         0     -1.25     -0.25         0  -1750.00
\end{verbatim}

En la Figura \ref{fig:eje_03} se muestra la región de factibilidad y la función objetivo en rojo para un valor $Z = 1400$. El código para determinar esta figura en matlab es:

\begin{verbatim}
hold off;

[x1,x2] = meshgrid(-100:1:600, -100:1:600);
f1 = 2*x1 + 3*x2 - 1200;
f2 = 2*x1 + x2 -1000;
f3 = 4*x2 - 800;
f4 = x1;
f5 = x2;

Z = 3*x1 + 4*x2 - 1400;

contour(x1, x2, f1, [0,0], 'k');
hold on;
contour(x1, x2, f2, [0,0], 'k');
contour(x1, x2, f3, [0,0], 'k');
contour(x1, x2, f4, [0,0], 'k');
contour(x1, x2, f5, [0,0], 'k');
contour(x1, x2, Z, [0,0], 'r');
xlabel('x_1');
ylabel('x_2');

\end{verbatim}

\begin{figure}[!htp]
   \centering
   \centerline{\psfig{figure=Capitulo_5/figuras/eje_03.jpg,width=14.0cm}}
   \caption{Ejemplo 3, área de Factibilidad en verde y en rojo función objetivo con $Z = 1400$.}
   \label{fig:eje_03}
\end{figure}
