Algoritmo PSO para la optimización de un controlador PI

En este Post trataremos una asunto que cada vez está cobrando más fuerza en el mundo de la Ingeniería, las técnicas inteligentes del aprendizaje. En concreto, os explicaré como desarrollar un algoritmo de optimización que sólo cambiando una función se puede aplicar a un sinfín de problemas. Dicho algoritmo es el PSO o Particle Swarm Optimization. En este caso el algoritmo será aplicado a la optimización de las constantes P (proporcional) e I (Integral) de un controlador tipo PI.

El algoritmo PSO es un método computacional que trata de optimizar un problema mejorando las soluciones propuestas, partículas, mediante la iteración de las mismas. Para discernir entre una buena y una mala solución se utiliza la función de coste, que consiste en asignar un valor a cada una de la soluciones con el fin de clasificarlas.  La función de coste será diferente para cada uno de los problemas a tratar. La solución con menor coste será considerada la mejor solución encontrada.

Para que el algoritmo funcione en primer lugar, es necesario la inicialización de todas las soluciones candidatas, en este caso se hará de forma aleatoria.  Después se calculan todos los costes,  y se inicializan los costes óptimos.  Una vez finalizada la inicialización se comienzan a iterar las soluciones, para ello, se calculan los costes de todas las partículas de la iteración y se actualizan los costes óptimos. Una vez finalizado el cálculo de los costes, se pasa a la actualización de la velocidad de cada una de las partículas para poder después actualizar las propias partículas. Estas actualizaciones se hacen en base a las siguientes expresiones:

Cabe destacar que los parámetros v', phi1 y phi2 se actualizarán en cada iteración de forma aleatoria entre unos máximos y mínimos fijados, mientras que la I (inercia) irá disminuyendo de forma proporcional al número de iteraciones que lleve la ejecución del algoritmo. Delta_t en cambio, es un parámetro constante a lo largo de toda la ejecución que marcará cuanto se alteran las partículas en cada iteración. 

  
clc; clear all;
//Algoritmo PSO
//Paso 1: Inicialización de las particulas: posicion inicial 
    //Calcular los costes para cada particula.
    //Calcular la posicion de cada particula.
    //Guardar el coste optimo de cada una de ellas.
    //Calcular la posicion optima de odo el enjambre.
    //Guardar el coste óptimo global.
n_itermax=100; 
n_particulas=60; 
n_var=2; 
Delta_t=1;
Inerciamax=5;
Inerciamin=0.1; 
Inercia = Inerciamax; 
phi_1_max=2; 
phi_2_max=2;
//Inicialización de las particulas:
x = random('uniform',0,150,1,n_particulas);
y = random('uniform',0,1,1,n_particulas);
x = [x; y];
x_optimos=x;
costes=ones(1,n_particulas);
for i=1:n_particulas
costes(i)=calcular_coste_controlador1(x(:,i));
end
costes_optimos=costes;
[coste_optimo_global,ref]=max(costes_optimos);
x_optimo_global=x_optimos(:,ref);

for n_iter=1:n_itermax
        for i=1:n_particulas
            //Calculo de coste de cada particula
           costes(i)=calcular_coste_controlador1(x(:,i));
            //Actualización del óptimo local
           if (costes(i)<=costes_optimos(i))
               costes_optimos(i)=costes(i);
               x_optimos(:,i)=x(:,i);
            //Actualización del óptimo global
               if costes_optimos(i)<=coste_optimo_global
                  coste_optimo_global=costes_optimos(i);
                  x_optimo_global=x_optimos(:,i);
               end
           end
        end

//Paso 2: Actualización de la velocidad de cada una de las particulas.
    //V_i(k)=I_i*V_i(k-1)+phi_u*(k_u_optimo(k)-X_i(k))+
    //phi_2*(X_global(k)-x_i(k))
    //Actualización aleatoria de los paramaetros para el 
    //calculo de las nuevas particulas
        phi_1=random('uniform',0,phi_1_max,1,1); 
        phi_2=random('uniform',0,phi_2_max,1,1);
        v=random('uniform',-0.1,0.1,n_var,n_particulas);
        for i=1:n_particulas
        //Calculo d ela velocidad
            v(:,i)=Inercia*v(:,i)+...
                phi_1*(x_optimos(:,i)-x(:,i))...
                +phi_2*(x_optimo_global-x(:,i));
            //Actualización de las posiciones:
            x(:,i)=x(:,i)+v(:,i)*Delta_t;
        end
//Extracción de constantes del controlador
Kp_optimo = x_optimo_global (1)
Ti_optimo = x_optimo_global (2)
//Creación de la función de trasnferencia del controlador
R = tf (Kp_optimo*[Ti_optimo,1],[Ti_optimo 0]);
//Creación de la funcion de transferencia a evaluar
H = tf (1,[1 ,1 ,1]);
//Creacion del lazo de control
Hla = R*H;
Hlc = feedback (Hla, 1);
//Parámetros de la simulación
tsimul = 20;
t = [0:0.001:tsimul]; //Vector de tiempo
yc = ones (1,length(t));//Vector de entrada
//Simualación y visualización del mejor resultado de esta iteración
y = lsim(Hlc,yc,t);
plot(t,y)
hold on
grid on
xlabel('Tiempo')
ylabel('Respuesta escalon')
100*n_iter/n_itermax %Visualización del timepo de ejecución
//Actualizacion de la inercia
Inercia = ((n_itermax-n_iter)/n_itermax)...
    *(Inerciamax-Inerciamin)+Inerciamin;
pause(0.1)
end
disp(x_optimo_global)  

A la hora de desarrollar este tipo de algoritmos la parte más difícil suele ser el diseño de una buena función de coste. No tiene un método fijo, por lo que es más bien imaginativa, eso sí, siempre se debe intentar que sea sencilla, que refleje bien problema y que tenga el mínimo coste computacional posible. 

Para este caso la función de coste planteada trata de igualar la FT analizada a una FT de referencia mediante el cálculo del error medio entre ambas. Además, se le han añadido otros parámetros a la función, como son el tiempo de establecimiento y el rebose máximo para conseguir llegar a la solución más rápidamente dotando de un coste más elevado a aquellas FT que tengan reboses y tiempos de establecimiento muy altos. Por lo tanto, la función de coste quedará de la siguiente manera:

Además de la función de coste en sí, se han implementado una serie de restricciones que harán que el coste de la FT evaluada sea infinito siempre y cuando no se cumplan. 

  
function [ coste ] = calcular_coste_controlador1( x )
//Extración de las constantes del controlador
Kp = x(1); 
Ti = x(2);
//Creacion de la funciond e trasferencia del PI
R = tf (Kp*[Ti,1],[Ti 0]);
//Creación de la funcion de transferencia a evaluar
H = tf (1,[1 ,1 ,1]);
//Creacion del lazo de control
Hla = R*H;
Hlc = feedback (Hla, 1);
//Creacion de funcion de transferencia de referencia
Href = tf (1,[1 1]);
//Parámetros de la simulación
tsimul = 20;
t = 0:0.01:tsimul; %Creación del vectro Tiempo
yc = ones (1,length(t)); %Creación vector de entrada
//Simaulación de la FT a analizar y la FT de referencia
yref = lsim (Href, yc, t);
y = lsim(Hlc,yc,t);
//Extracción de datos de la FT a analizar
datos = stepinfo(Hlc,tsimul,1);
if Kp < 0 || Ti < 0  %Restricción 1
    coste = inf;
else
    if datos.PeakTime > 5 %Restricción 2
        coste = inf; 
    else 
        if datos.Overshoot > 60
            coste = inf;
        else
            //Función de coste
            coste = mean(yref-y)/1 + (datos.Overshoot)/60 + (datos.SettlingTime)/15;
        end  
    end
end
end

Después de varias ejecuciones con diferentes ajustes en los parámetros del algoritmo estos son los resultados obtenidos con 100 iteraciones y 60 partículas.

Analizando los resultados se puede concluir que se ha llegado a un resultado final bastante satisfactorio donde el tiempo de establecimiento ronda los 8 segundos y donde apenas hay un 1% de rebose.