OPTIMIZACION Y

CONTROL
DO IT YOURSELF
AYOUB EL MKADEM BOUKLOUZ
EJERCICIO 1

Antes de resolver un ejercicio numérico un ingeniero debe de tener la capacidad critica de

sentarse y estudiar la función antes de hacer nada, es decir primero hay que jugar con la
función antes de querer resolverla, ver los limites estudiarla, sentir como es…

Esto mencionado es lo que diferencia un buen analista de otro normal, y esto se ve reflejado
en la resolución de la función del ejercicio 1. Todos sabemos darle a un botón y resolver un
script o copiar cuatro líneas de código y cambiar la función pero estudiar y analizar los
resultado son cosas diferentes y es lo que a la larga va a hacer que una persona destaque
sobre otra, LA CAPACIDAD CRITICA Y ANALISTA DE UNA PERSONA, empecemos:

En esta función si nos ponemos a lo loco con el Matlab a resolverlo con las dos funciones que
la profesora ha proporcionado, fminsearch y fminunc, va a ver lo siguiente:

Para el primer punto de iteración [ 1 1]

Para fminunc:
Para fminsearch

Donde una representación de los valores para cada iteración viene dada por la siguiente
gráfica, donde se ve como la función converge a un determinado valor:
Si hacemos una representación de la función dada cerca del punto local donde están
encontrando mis funciones el mínimo, vemos lo siguiente:
Vemos que efectivamente se trata de un mínimo, pero no sabemos si realmente es el mínimo
local o un mínimo global y eso lo sabremos si ampliamos un poco más nuestro espectro de
representación, pero aun más importante si analizamos la función como decíamos en el primer
párrafo antes de hacer cualquier cosa, para tener una idea de qué es lo que estamos haciendo
para saber por donde irán los tiros y no solamente fiarse de cualquier resultado, lo que
menciono es super importante para cálculos numéricos en solidos o fluidos, uno no puede
fiarse de cualquier cosa debe de analizar los resultados.

La función dada viene dada por:

− cos(𝑥) · cos(𝑦)
𝑓(𝑥, 𝑦) = 2 +(𝑦−𝜋)2
e(x1−π)

Si nos fijamos en la función, vemos que existe un mínimo global y es para los valores de

𝑥=𝜋
𝑦=𝜋
𝑓(𝜋, 𝜋) = −1
Ese es el mínimo valor que puede alcanzar esa función, porque tenemos que el numerador
será periódico, y el denominador si crece mucho el valor se me irá a cero por la exponencial,
pero si los valores de las incógnitas son igual a pi, tenemos que el denominador es 1, nunca
puede ser negativo. Como estamos minimizando y no maximizando debemos buscar un valor
de los cosenos que me maximice la función y se mantenga negativa, es decir que ambos
cosenos valgan el máximo valor que será 1 o -1. Pero los cosenos valen 1 cuando x e y son 0,
pero con esos valores tendríamos que el denominador es mayor de uno, por lo tanto, el valor
de la función tiende a 0, Entonces el valor mínimo de la función se alcanza cuando la x y la y
son pi y por lo tanto el valor mínimo es -1. Para verificar este análisis vamos a ampliar el
espectro de la supercicie representada y vamos a ver si existe un mínimo global.
Y efectivamente se cumple que el punto que hemos encontrado antes es un monimo local que
se encontraba en una zona plana, donde se alacanza que el gradiente es cero y nos devuelve el
mínimo que encontramos, para los siguientes puntos vemos que la función no hace nada y nos
devuelve los mismos puntos de iteración, eso quiere decir que se trtaa de una zona plana en la
que el gradiente es cero. Por lo tanto y como dijimos tendremos un mínimo global que alcanza
el valor de -1 para los valores de pi tanto para los valores de x como de y.

Uno puede pensar que se trata de la tolerancia, pero hemo intentado cambiar la tolerancia
que viene por defecto, pero no cambia nada, sigue encontrando para el punto inicial [1 1] el
mismo valor, en la función fminunc, nos da un aviso de que se trata de un mínimo local,
incluso si cambiamos las tolerancias los resultados no cambian ni para fminsearch ni para
fminunc, por lo tanto, el valor de 1.305 para la x y la es el mínimo local que encuentra para un
punto inicial de [1 1].

Con respecto a los demás puntos iniciales:

V02 = [ 20 15]

fminunc fminsearch
Efectivamente como hemos comentado no se mueve del punto de iteración porque está en
una zona plana.

V03 = [-2 5]

Fminunc fminsearch

En este caso el fminunc no se mueve por lo que comentamos del gradiente y el fminsearch
encuentra una mejor solución y se mueve a otro punto, aunque sin llegar a la solución del
mínimo global del problema.

Los tiempos que tarda cada uno son muy diferente, para este ultimo caso, el tiempo que ha
tardado en ejecutarse el fminunc es de 0.08 segundo y el fminsearch es de 1.3 segundos.
Vemos que la función fminunc es unas diez veces más rápida o más.

Además, hemos usado la función de islocalmin para ver los mínimos que presenta la función,
se usa mucho en señales, podemos encontrar las posiciones de los elementos que son
mínimos locales, si aumentamos el espectro vemos que encuentra varios mínimos locales
ahora representamos cuales son y los valores de la función en esos puntos
Como vemos para x1 y x2 igual a pi o cerca tenemos que la función vale -0.998 que es el valor
mínimo de los representados y por lo tanto es un mínimo global. Y los demás puntos también
encontramos el pinto detectado para el punto [1 1] que encuentra el mínimo local para el
valor de [1.3 1.3]. Por lo tanto, tenemos muchos mínimos locales y un mínimo global. Es
verdad que uno podría no considerar estos puntos como mínimos locales porque no se trata
de una función que una vez representada en 3D se vean a ojo como se ve en otros ejercicios,
pero para mí personalmente la representación de la superficie cerca del primer punto,
muestra claramente que en esa región hay un mínimo local, además si empiezo en un punto y
el fminsearch consigue encontrar una zona donde la función vale menos pues es muestra de
que esa zona no es del todo plana donde estamos buscando, ¿o se tratará de errores de
aproximación por truncamiento?, si analizamos la función vemos que si cambiamos los valores
de x1 y x2 tendremos valores diferentes, muy pequeños pero diferentes, si truncásemos esto,
pues si que tendríamos una zona plana pero para mi no es del todo plana sino que tiene
llamémoslo ‘’baches’’.

Hemos probado con otros puntos iniciales y uno de ellos es un punto cercano al mínimo global.
Los resultados obtenidos son los siguiente:
Fminunc fminsearch

Efectivamente ambos métodos encuentran el mínimo global en el punto [pi,pi].

fminunc tarda mucho menos comparar los tiempos y el número de iteraciones también es
menor

Por ello vemos que la información del gradiente ayuda a alcanzar la solución más rápido

En el script de Python puede encontrar los miles de posibilidades que hemos verificado, con
diferentes tolerancias, diferentes algoritmos de cálculo, hasta hemos comparado entre darle al
programa un gradiente analítico, pero no hemos notado una mejoría al darle un gradiente
comparado con el gradiente numérico que calcula con el método de bfgs, no explicamos cada
detalle en este informe porque sería de 100 páginas.

De hecho, en los siguientes problemas seré breve y contestaré solamente a las preguntas que
se plantean. Y en los códigos muchas cosas que hemos probado las dejamos comentadas para
no calcular todo a la vez pero que puede verificar los resultados quitando el ‘%’ de esas líneas.
EJERCICIO 2

Probamos con mu igual a cero todo y lo que encuentra es el mínimo [6 7] que sería como
resolver el problema sin restricciones. La representación 3D sería:
Se va a resolver el problema como una función única con penalty, es decir a la función objetivo
le añadimos las restricciones multiplicados por una constante mu. Las restricciones tenemos
que ponerlas de la forma 𝑔(𝑥) ≤ 0. Para el método cuadrático tenemos el siguiente gráfico si
establecemos cero para los cuatro mu de las cuatro restricciones:

La región donde tenemos que encontrar la solución es dentro del rombo que forman las cuatro
líneas que se corresponden con las cuatro restricciones. Vemos que para un valor de mu igual
a cero para todos, vemos que la solución tiende al mínimo global del problema, porque si
usamos mu igual a cero es como si no usáramos las restricciones.

A continuación, empezamos a subir hasta los valores, de mu, y nos damos cuenta que la
frontera que está activa es la correspondiente a la tercera frontera:

𝑥1 + 𝑥2 − 7 ≤ 0
Por lo tanto, el valor de mu3 del código es el que hace que la solución se quede en la frontera.

La restricción 3 está activa, es la que manda en la solución y es con la que tenemos que lidiar
para que cumpla la función las demás se cumplen perfectamente.

A continuación, presentamos algunas imágenes para diferentes valores de mu:

mu1 = 100 ; mu2 = 100; mu3 = 100; mu4 = 100;

LA solución óptima que encuentra dentro de la región es:

Para mu1 = 100 ; mu2 = 100; mu3 = 1; mu4 = 100;

Vemos que no está dentro de la región, por lo tanto, no es válida.

Después de muchas iteraciones, hemos visto que el valor de mu3 correspondiente a la

restricción 3 es la que rige que la solución esté dentro de la región factible del problema.
Para mu1 = 1; mu2 = 1; mu3 = 100; mu4 = 1;

En este anterior vemos que aunque los demás ‘mu’ sean pequeños, con que sea grande el
mu3, hace que la solución esté dentro de la región factible.

Para los métodos de barrera hemos programado tanto el de barrera inversa como el
logarítmico. Para estos ha sido muy difícil encontrar un mu que cumpla con las condiciones y
dé una solución real. Por lo tanto, para el método de barrera inversa hemos creado una
estructura de cuatro bucles para encontrar los cuatro mu que permita que la solución esté en
la región factible.

Dentro de la estructura de los bucles hemos puesto una condición que cuando resuelva el
problema de optimización verifique que esos puntos están dentro de la región factible y
además que el valor de la función no sea igual a infinito o menos infinito. Los puntos que
cumplen con la región factible y que su función sea finita, los guardamos en diferentes
vectores. Además, el initial guess, hemos escogido un punto que esté dentro de la región, en
este caso el [2 2] porque los métodos de barrera son necesario que el punto inicial esté dentro
de la región para que no salga de la barrera.
Para los métodos de barrera no hemos tenido mucho éxito la verdad, están programados en el
mismo script, hemos probado para mucho mu pero no hemos conseguido grandes resultados.

Para este método hemos obtenido los siguientes valores mínimos:

Sabemos que no es la mejor solución, pero para ello tendríamos que aumentar la resolución
de las constantes mu, hemos probado pero el problema tarda mucho lo que no es
recomendable, por lo tanto nos quedaríamos con esta solución que está dentro de la región
factible pero no es la mejor.

Con respecto al método de barrera logarítmico, lo programamos con la misma estructura que
el inverso y nos da:
Un resultado parecido, pero que sin duda el valor optimo fue alcanzado con el método de
penalty cuadrático, ya que es más robusto y tampoco necesitamos darle el punto dentro de la
región factible, que en un problema muy complejo es posible que no la conozcamos. No
hemos comparado los tiempos porque hemos creado un bucle muy largo para los métodos de
barrera por lo tanto no tendría sentido hablar de tiempos en este apartado, pero claramente el
cuadrático es más eficiente y más rápido.

EJERCICIO 3

Para resolver el problema con como problema sin restricciones, el problema lo vamos a
resolver con el método de penalty cuadrático. Para ello tenemos cuatro restricciones:

𝑥1 + 𝑥2 − 5 ≤ 0
−𝑥1 ≤ 0
𝑥1 ≤ 3
−𝑥2 ≤ 0

La forma de la superficie es la que se presenta a continuación:

En el siguiente ejercicio programamos dos scripts, uno que lo resolvemos con los métodos de
sqp, que es mediante la función fminunc y la otra con métodos de penalty cuadráticos. En
ambos scripts conseguimos un resultado de:

SQP Penalty

Vemos que ambos métodos llegan al mismo valor y que coincide con la zona que presentamos
en la superficie 3D donde se alcanza el mínimo, de [3 2] donde la función vale -31. Pero la gran
diferencia está en la velocidad de cálculo. Claramente el método SQP es 10 veces más rápido
que los métodos de penalty y esto en un problema muy complejo donde el tiempo es
importante claramente el ganador sería el método SQP.
EJERCICIO 4

La superficie en 3D tiene la siguiente forma:

Claramente tiene muchos mínimos locales, pero para detectar el absoluto usaremos tres
métodos, el multistart, el hibrido y el differential evolution.

Con la función multistart vemos que la función encuentra todos los mínimos locales y escoge el
menor de ellos. Después de las 200 evaluaciones nos sale la mejor.
Encuentra la mejor solución para los puntos aproximadamente [ 0 0 ], que si sustituimos en la
función original sale un valor de cero. A continuación, presentamos un histograma que recoge
las selecciones de los valores que ha ido escogiendo y finalmente un contorno donde los
puntos rojos es donde encuentra el mínimo local.
Como vemos y hemos comentado al principio, hay muchos mínimos locales y el programa
selecciona el mejor de ellos que en este caso es el [0 0]. Pero el programa tarda como 6
segundos en encontrar el mínimo global.

Si aumentamos el espectro [ -50 50]:

Lo que pasa aquí es que, al tener un espectro tan grande, para 200 evaluaciones, pierde mucha
precisión con respecto al anterior que era un espectro más pequeño, aquí tendríamos que
añadir el numero de evaluaciones si queremos que tenga valores más precisos, pero aun así se
obtienen valores razonables.

Con el método del DE, obtenemos los siguientes resultados:

Vemos que encuentra un mínimo global en un punto cercano al [0 0 ] y que el valor es

aproximadamente 0 de la función, aquí podríamos jugar con las estrategias seleccionadas,
hemos seleccionado la 7 y también podríamos cambiar parámetros como la F o CR o el numero
de iteraciones máxima, sería probar, pero como obtenemos un punto cercano al valor del
mínimo global lo damos por válido, además el tiempo que ha tardado es de 0.25 segundos, es
decir 25 veces menos que con el método del multistart.

Probamos otro método que es el híbrido y obtenemos:

Con el método hibrido hemos obtenido resultados similares, un mínimo global cerca del [0 0+
y de un valor de 0 para la función pero que es mucho más rápido que el multistart y que el DE.
Una 100 veces más rápido que el multistart y unas 5 más que el DE.
EJERCICIO 5

Este problema los hemos resuelto con tres métodos para encontrar el Pareto óptimo de estas
dos funciones para las restricciones que nos dan.

Método de pesos:

Como este problema intenta resolver una única función dándole una serie de pesos a cada
una, además hemos añadido las restricciones de forma que resolvemos el optimo con métodos
de penalización cuadráticos. Es decir, las restricciones son:

−𝑥1 ≤ 0
𝑥1 − 5 ≤ 0
−𝑥2 ≤ 0
𝑥2 − 5 ≤ 0

Nos queda en el Matlab de la forma:

Resolvemos para una serie de pesos y nos vamos quedando con la solución optima para cada
peso dado, que posteriormente representarán la frontera del óptimo de Pareto.

Una vez resuelto el script nos queda:

Los puntos azules representan el óptimo de Pareto y cualquiera de ellos sería una solución
buena a mi problema,

A continuación, presentamos una tabla con los resultados de los valores x1, x2, f1 y f2
correspondientemente:
Como vemos la x1 y x2 van de 0 a 5, eso quiere decir que se han considerado las restricciones
añadidas y el valor de f1 y f2 tienen sentido ya que para el máximo de f1 se alcanza el mínimo
de f2 y viceversa. Un punto óptimo podría ser el [ 1.1702 1.1702] que corresponde a los
valores de f1 y f2 [10.95 29.33].

Ahora resolvemos el mismo problema con dos métodos diferentes, el e-constrained y por
nsga2.

Método e-constrained:

Este método fija un valor de f2, establecemos un delta y para cada tramo vamos resolviendo el
óptimo. Este código y el anterior han sido programados desde cero por mí, no se ha copiado
de ningún sitio, se han programado siguiendo las lecciones de clase, por ello no explicaremos
en detalle que hace.

Ejecutamos el programa y obtenemos la siguiente gráfica que coincide exactamente con la

obtenida en el método de pesos.

También se han tenido en cuenta las restricciones como en el caso anterior al resolver el
problema de optimización una vez se fija f2. Algunos resultados correspondientes a x1, x2, f1 y
f2 del frente óptimo de Pareto es:
Si cogemos un punto próximo al anterior obtenemos el [ 1.127 1.1269] que se obtiene
resultados de f1 y f2 de [11.68 28.75] que son bastante parecidos a los anteriores. Por lo tanto,
concluimos que este método ha funcionado perfectamente.

Método NSGAII:

En este script hemos adaptado nuestro problema al script que nos ha proporcionado la
profesora, por lo tanto, no era muy complejo al resolver, para comparar vemos que el frente
de Pareto sale parecido a los anteriores y además hemos buscado esos puntos que hemos
comparado anteriormente y nos da un valor idéntico, por lo tanto hemos obtenido resultados
bastante buenos para los tres métodos. En el Matlab se puede ver el trabajo de cada script y el
último como lo hemos adaptado para nuestro problema.

Efectivamente la f1 = 11.1 y f2 = 29.31 se parece mucho a los puntos seleccionados de los

frentes de Pareto anteriores.

Para ejecutar este archivo, entrar en la carpeta del ejercicio 5, NGSA2 y buscar la función
doit_5_NSGA2.m

Para ellos tiempos de cálculo, tenemos que el e-constrained y el método del peso tardan
aproximadamente lo mismo 0.6 segundos, en cambio este ultimo tarda un poco más porque
tenemos activada la representación. Pero este método es más preciso junto con el e
constrained porque el método de los pesos solo me obtiene el Pareto optimo en caso de que
sea convexa, en cambio estos dos últimos cualquiera.