TÓPICOS EN SISTEMAS EXPERTOS E INTELIGENCIA ARTIFICIAL 1. Introducción En el presente trabajo se desarrolló el problema 7.4 del libro “Power Generation, Operation and Control” [1] mediante el método metaheurístico de Optimización por Enjambre de Partículas. Para esto se implementó un programa en el software matemático MATLAB. Como primer paso, se desarrolló el ejemplo 7.E del libro [1], a fin de validar la coherencia de los resultados obtenidos mediante el modelo implementado en Matlab de este método metaheurístico. Cabe mencionar que en el libro de referencia se desarrolla el ejemplo 7.E mediante el método matemático de programación dinámica. 2. Ejemplo 7.E 2.1 Enunciado del Ejemplo 7.E Realizar el despacho económico para el sistema de la figura 1. Figura 1 Sistema del ejemplo 7.E.  Características de la planta de vapor “S”: Figura 2 Función de costos incrementales.  Características de la planta hidroeléctrica “H”: Figura 3 Función de caudal vs Potencia generada. TÓPICOS EN SISTEMAS EXPERTOS E INTELIGENCIA ARTIFICIAL Volumen inicial=Volumen final=10000 acre-ft 6000≤Volumen≤18000 Caudal natural=1000 acre-ft/h  Característica de la demanda Tabla 1 Característica de la demanda 2.2 Formulación Matemática del problema 𝑗=6 𝑀𝑖𝑛 { ∑ 𝐹(𝑃𝑠(𝑗) ) × 4} 𝑗=1 s.a. 𝑃𝑠(𝑗) + 𝑃𝐻(𝑗) = 𝐷𝑗 200 ≤ 𝑃𝑠(𝑗) ≤ 1200 0 ≤ 𝑃𝐻(𝑗) ≤ 200 6000 ≤ 𝑉(𝑗) ≤ 18000 𝑉(0) = 10000 𝑉(6) = 10000 𝑉(𝑗) = 𝑉(𝑗−1) + 4(1000 − 𝑞(𝑗) ) 𝑞(𝑗) = 260 + 10 × 𝑃ℎ(𝑗) 3. Modelamiento en Matlab Para simplificar el código MATLAB implementamos dos funciones “ofun(x)” y “Resultados(x)”. La primera función calcula el valor de la función objetivo para el vector “x” ingresado. La segunda función grafica los resultados de generación y volúmenes para el vector “x” ingresado. Las restricciones se modelaron con una penalización de R/. 10 000, con el fin de que estos resultados tengan un valor de función objetivo muy alto respecto al resultado óptimo. Esto con el fin de evitar que sean consideradas como mejores posiciones durante el desarrollo del método heurístico PSO. El código implementado en Matlab se presenta en el anexo 1. TÓPICOS EN SISTEMAS EXPERTOS E INTELIGENCIA ARTIFICIAL 3.1 Análisis de los resultados obtenidos En las figuras 4 y tabla 2 se presentan los resultados del ejemplo 7.E obtenidos mediante el método matemático de Programación Dinámica. En la figura 4 se presenta el volumen del embalse de la planta hidroeléctrica asociada a la operación óptima del sistema. Mientras en la tabla 2 se resalta el valor óptimo de la función objetivo. Figura 4 Comportamiento del volúmen del embalse de la planta hidroeléctrica asociado a la operación óptima del sistema. Tabla 2 Función Objetivo optima. A continuación se presenta los resultados obtenidos mediante el método heurístico de Enjambre de Partículas para el mismo problema. En la figura 5, se presenta la potencia generada óptima por cada planta y para cada período. En la figura 6, se presenta el comportamiento del volumen del embalse de la planta hidroeléctrica asociada a la operación óptima del sistema. En la tabla 3, presentamos la función objetivo obtenido con el método PSO. Nótese que los resultados obtenidos del método heurístico son coherentes con el resultado obtenido del método matemático. Además, es importante hacer notar que los resultados son diferentes en 0.14 % esto debido a que no necesariamente un método heurístico llega a la solución óptima global del problema, sin embargo es importante precisar que para este caso se obtuvo un resultado bastante aproximado al optimo real. Programación Dinámica PSO Diferencia de resultados (%) Función Objetivo Óptima 81738 81853 0.14% Tabla 3 Funciones Objetivos optimas obtenidas con metodo matemático ymetodo heurístico. TÓPICOS EN SISTEMAS EXPERTOS E INTELIGENCIA ARTIFICIAL Resumen de Resultados - Generación 1000 Planta de Vapor (MW ) Planta Hidroeléctrica (MW ) 900 800 700 (MW) 600 500 400 300 200 100 0 1 2 3 4 Periodo "j" 5 6 Figura 5 Generación de la planta hidroeléctrica y Generación de la planta de vapor para cada período. 1.5 x 10 4 Resumen de Resultados - Volumen(j) 1.4 Volúmen (acre-ft) 1.3 1.2 1.1 1 0.9 0.8 0.7 0.6 1 2 3 4 Periodo "j" 5 6 7 Figura 6 Volúmen del embalse de la planta hidroeléctrica para cada período. Luego, podemos decir que se obtuvo resultados coherentes de nuestro modelo implementado en Matlab, basados en la verificación realizada con los resultados del ejemplo 7.E. Por lo tanto, luego de haber verificado la validez de nuestro modelo, se procedió a solucionar el problema 7.4. TÓPICOS EN SISTEMAS EXPERTOS E INTELIGENCIA ARTIFICIAL 4. Problema 7.4 4.1 Enunciado del Problema 7.4 TÓPICOS EN SISTEMAS EXPERTOS E INTELIGENCIA ARTIFICIAL 4.2 Formulación Matemática del problema 𝑗=6 𝑀𝑖𝑛 { ∑ 𝐹(𝑃𝑠(𝑗) ) × 4} 𝑗=1 s.a. 𝑃𝑠(𝑗) + 𝑃𝐻(𝑗) = 𝐷𝑗 200 ≤ 𝑃𝑠(𝑗) ≤ 1200 0 ≤ 𝑃𝐻(𝑗) ≤ 2000 (0.9 + 𝑉(𝑗) + 𝑉(𝑗−1) ) 200000 6000 ≤ 𝑉(𝑗) ≤ 18000 𝑉(0) = 10000 𝑉(6) = 10000 𝑉(𝑗) = 𝑉(𝑗−1) + 4 {1000 − (260 + 10𝑃𝐻(𝑗) (1.1 − 𝑉(𝑗) + 𝑉(𝑗−1) ))} 200000 4.3 Modelamiento en Matlab Se actualizó los nuevos datos en el modelo implementado en 2.3. El código implementado en Matlab se presenta en el anexo 2. 4.4 Análisis de los resultados obtenidos A continuación se presenta los resultados obtenidos mediante el método heurístico de Enjambre de Partículas. En la figura 7, se presenta la potencia generada óptima por cada planta y para cada período. En la figura 8, se presenta el comportamiento del volumen del embalse de la planta hidroeléctrica asociada a la operación óptima del sistema. En la tabla 4, presentamos la función objetivo obtenido con el método PSO. Además se presenta la generación óptima de las plantas para cada período “j”. Es importante hacer notar que no necesariamente un método heurístico llega a la solución óptima global del problema, sin embargo es importante precisar que en el punto 3.4 se obtuvo un resultado bastante aproximado al óptimo real. Período 1 2 3 4 Ps 600 1000 545.7 500 PH 0 0 354.3 0 Costo total de Operación (Función Objetivo Óptima)= R/. 89 401 5 200 200 Tabla 4. Generación Óptima de las planta por cada período “j” y Función Objetivo Optima. 6 300 0 TÓPICOS EN SISTEMAS EXPERTOS E INTELIGENCIA ARTIFICIAL Resumen de Resultados - Generación 1000 Planta de Vapor (MW ) Planta Hidroeléctrica (MW ) 900 800 700 (MW) 600 500 400 300 200 100 0 1 2 3 4 Periodo "j" 5 6 Figura 7 Generación óptima de la planta de vapor y de la planta hidroeléctrica. 1.8 x 10 4 Resumen de Resultados - Volumen(j) 1.6 Volúmen (acre-ft) 1.4 1.2 1 0.8 0.6 0.4 1 2 3 4 Periodo "j" 5 6 7 Figura 8 Volúmen del embalse de la planta hidroeléctrica asociada a la operación óptima del sistema para cada período “j”. Suárez Limache José Alberto TÓPICOS EN SISTEMAS EXPERTOS E INTELIGENCIA ARTIFICIAL Anexo 1 Código MATLAB – Ejemplo 7.E TÓPICOS EN SISTEMAS EXPERTOS E INTELIGENCIA ARTIFICIAL v(i,j)=w*v(i,j)+c1*rand()*(pbest(i,j)x(i,j))+c2*rand()*(gbest(1,j)-x(i,j)); end end %Actualización de la posición de las "n" particulas. for i=1:n for j=1:m x(i,j)=x(i,j)+v(i,j); end end %Corrección de las posiciones de las "n" particulas en caso violen algún límite. for i=1:n for j=1:m if x(i,j)<LB(j) x(i,j)=LB(j); elseif x(i,j)>UB(j) x(i,j)=UB(j); end end end %Calculamos la función objetivo para las nuevas posiciones de las "n" partículas. for i=1:n f(i,1)=ofun(x(i,:)); end %Se grabará las mejores posiciones por cada partícula. for i=1:n if f(i,1)<f0(i,1) pbest(i,:)=x(i,:); %pbest: mejores posiciones. f0(i,1)=f(i,1); %f0: función objetivo de las mejores posiciones. end end [fmin,index]=min(f0); %Mejor posición de las mejores posiciones locales %Actualizando la mejor posición global if fmin<fmin0 gbest=pbest(index,:); fmin0=fmin; end %Presentamos los resultados de la mejor partícula. if ite==1 disp('Iteración Mejor partícula Función Objetivo'); end disp(sprintf('%8g %8g %20g',ite,index,fmin0)); ite=ite+1; end %--------------------------------------------------------------------Funcion_Objetivo_Optima(run)=ofun(gbest); Posicion_Optima(run,:)=gbest; disp('----------------------------------------') end % pso main program-----------------------------------------------------end TÓPICOS EN SISTEMAS EXPERTOS E INTELIGENCIA ARTIFICIAL disp('**************************************************************** **'); disp('Resumen de resultados'); [Funcion_Objetivo_Optima_Final,bestrun]=min(Funcion_Objetivo_Optima); Funcion_Objetivo_Optima_Final bestrun Variables_optimas=Posicion_Optima(bestrun,:) disp('**************************************************************** **'); toc %Graficando Volumen del embalse Resultados(Variables_optimas); TÓPICOS EN SISTEMAS EXPERTOS E INTELIGENCIA ARTIFICIAL function f=ofun(x) %Esta función calculará la función objetivo de la posición "x" %Función Obejtivo "Min" of=(700+(4.8*x(1))+(x(1)*x(1)/2000))*4+(700+(4.8*x(2))+(x(2)*x(2)/2000 ))*4+(700+(4.8*x(3))+(x(3)*x(3)/2000))*4+(700+(4.8*x(4))+(x(4)*x(4)/20 00))*4+(700+(4.8*x(5))+(x(5)*x(5)/2000))*4+(700+(4.8*x(6))+(x(6)*x(6)/ 2000))*4; %Calculamos la generación de la planta hidroeléctrica. Ph1=600-x(1); Ph2=1000-x(2); Ph3=900-x(3); Ph4=500-x(4); Ph5=400-x(5); Ph6=300-x(6); V0=10000; %Calculamos los volúmenes. if Ph1==0 V1=V0+4*(1000); %q1=0; else V1=V0+4*(1000-(260+10*Ph1)); %q1=260+10*Ph1; end if Ph2==0 V2=V1+4*(1000); %q2=0; else V2=V1+4*(1000-(260+10*Ph2)); %q2=260+10*Ph2; end if Ph3==0 V3=V2+4*(1000); %q3=0; else V3=V2+4*(1000-(260+10*Ph3)); %q3=260+10*Ph3; end if Ph4==0 V4=V3+4*(1000); %q4=0; else V4=V3+4*(1000-(260+10*Ph4)); %q4=260+10*Ph4; end if Ph5==0 V5=V4+4*(1000); %q5=0; else V5=V4+4*(1000-(260+10*Ph5)); %q5=260+10*Ph5; end if Ph6==0 V6=V5+4*(1000); %q6=0; else V6=V5+4*(1000-(260+10*Ph6)); %q6=260+10*Ph6; end %Restricciones c0=[]; c0(1)=6000-V1; c0(2)=V1-18000; c0(3)=6000-V2; c0(4)=V2-18000; c0(5)=6000-V3; c0(6)=V3-18000; c0(7)=6000-V4; c0(8)=V4-18000; c0(9)=6000-V5; % % % % % % % % % <=0 <=0 <=0 <=0 <=0 <=0 <=0 <=0 <=0 TÓPICOS EN SISTEMAS EXPERTOS E INTELIGENCIA ARTIFICIAL c0(10)=V5-18000;% <=0 c0(11)=9999-V6; % <=0 c0(12)=V6-10001;% <=0 %Creamos la matriz "c" con "ceros" en caso cumpla las restricciones y "unos" en caso viole las restricciones. for i=1:length(c0) if c0(i)>0 c(i)=1; else c(i)=0; end end %Penalizamos las violaciones a las restricciones. penalicacion=100000; %Función Objetivo de la posición "x" f=of+penalicacion*sum(c); TÓPICOS EN SISTEMAS EXPERTOS E INTELIGENCIA ARTIFICIAL function f=Resultados(x) %Graficando reusltados de generación D=[600 1000 900 500 400 300]; %Actualizar la demanda for i=1:length(x) gen(i,1)=x(i); end for e=1:length(x) gen(e,2)=D(e)-x(e); end bar(gen,'stack') title('Resumen de Resultados - Generación') xlabel('Periodo "j"') ylabel('(MW)') legend('Planta de Vapor (MW)','Planta Hidroeléctrica (MW)') grid on %Graficando Volúmenes V(1)=10000; %Calculamos los volúmenes. for i=1:length(x) if gen(i,2)==0 V(i+1)=V(i)+4*(1000); %q1=0; else V(i+1)=V(i)+4*(1000-(260+10*gen(i,2))); %q1=260+10*Ph1; end end figure plot(V) title('Resumen de Resultados - Volumen(j)') xlabel('Periodo "j"') ylabel('Volúmen (acre-ft)') grid on TÓPICOS EN SISTEMAS EXPERTOS E INTELIGENCIA ARTIFICIAL Anexo 2 Código MATLAB – Problema 7.4 TÓPICOS EN SISTEMAS EXPERTOS E INTELIGENCIA ARTIFICIAL v(i,j)=w*v(i,j)+c1*rand()*(pbest(i,j)x(i,j))+c2*rand()*(gbest(1,j)-x(i,j)); end end %Actualización de la posición de las "n" particulas. for i=1:n for j=1:m x(i,j)=x(i,j)+v(i,j); end end %Corrección de las posiciones de las "n" particulas en caso violen algún límite. for i=1:n for j=1:m if x(i,j)<LB(j) x(i,j)=LB(j); elseif x(i,j)>UB(j) x(i,j)=UB(j); end end end %Calculamos la función objetivo para las nuevas posiciones de las "n" partículas. for i=1:n f(i,1)=ofun(x(i,:)); end %Se grabará las mejores posiciones por cada partícula. for i=1:n if f(i,1)<f0(i,1) pbest(i,:)=x(i,:); %pbest: mejores posiciones. f0(i,1)=f(i,1); %f0: función objetivo de las mejores posiciones. end end [fmin,index]=min(f0); %Mejor posición de las mejores posiciones locales %Actualizando la mejor posición global if fmin<fmin0 gbest=pbest(index,:); fmin0=fmin; end %Presentamos los resultados de la mejor partícula. if ite==1 disp('Iteración Mejor partícula Función Objetivo'); end disp(sprintf('%8g %8g %20g',ite,index,fmin0)); ite=ite+1; end %--------------------------------------------------------------------Funcion_Objetivo_Optima(run)=ofun(gbest); Posicion_Optima(run,:)=gbest; disp('----------------------------------------') end % pso main program-----------------------------------------------------end TÓPICOS EN SISTEMAS EXPERTOS E INTELIGENCIA ARTIFICIAL disp('**************************************************************** **'); disp('Resumen de resultados'); [Funcion_Objetivo_Optima_Final,bestrun]=min(Funcion_Objetivo_Optima); Funcion_Objetivo_Optima_Final bestrun Variables_optimas=Posicion_Optima(bestrun,:) disp('**************************************************************** **'); toc %Graficando Volumen del embalse Resultados(Variables_optimas); TÓPICOS EN SISTEMAS EXPERTOS E INTELIGENCIA ARTIFICIAL function f=ofun(x) %Esta función calculará la función objetivo de la posición "x" %Función Obejtivo "Min" of=(770+(5.28*x(1))+(x(1)*x(1)*0.00055))*4+(770+(5.28*x(2))+(x(2)*x(2) *0.00055))*4+(770+(5.28*x(3))+(x(3)*x(3)*0.00055))*4+(770+(5.28*x(4))+ (x(4)*x(4)*0.00055))*4+(770+(5.28*x(5))+(x(5)*x(5)*0.00055))*4+(770+(5 .28*x(6))+(x(6)*x(6)*0.00055))*4; %Calculamos la generación de la planta hidroeléctrica. Ph1=600-x(1); Ph2=1000-x(2); Ph3=900-x(3); Ph4=500-x(4); Ph5=400-x(5); Ph6=300-x(6); V0=10000; %Calculamos los volúmenes. if Ph1==0 V1=V0+4*(1000); %q1=0; else V1=(5000*V0+14800000-220000*Ph1+Ph1*V0)/(5000-Ph1); %q1=260+10*Ph1*(1.1-(V0+V1)/200000); end if Ph2==0 V2=V1+4*(1000); %q2=0; else V2=(5000*V1+14800000-220000*Ph2+Ph2*V1)/(5000-Ph2); %q2=260+10*Ph2*(1.1-(V1+V2)/200000); end if Ph3==0 V3=V2+4*(1000); %q3=0; else V3=(5000*V2+14800000-220000*Ph3+Ph3*V2)/(5000-Ph3); %q3=260+10*Ph3*(1.1-(V2+V3)/200000); end if Ph4==0 V4=V3+4*(1000); %q4=0; else V4=(5000*V3+14800000-220000*Ph4+Ph4*V3)/(5000-Ph4); %q4=260+10*Ph4*(1.1-(V3+V4)/200000); end if Ph5==0 V5=V4+4*(1000); %q5=0; else V5=(5000*V4+14800000-220000*Ph5+Ph5*V4)/(5000-Ph5); %q5=260+10*Ph5*(1.1-(V4+V5)/200000); end if Ph6==0 V6=V5+4*(1000); %q6=0; else V6=(5000*V5+14800000-220000*Ph6+Ph6*V5)/(5000-Ph6); %q6=260+10*Ph6*(1.1-(V5+V6)/200000); end %Restricciones c0=[]; c0(1)=6000-V1; % <=0 TÓPICOS EN SISTEMAS EXPERTOS E INTELIGENCIA ARTIFICIAL c0(2)=V1-18000; % <=0 c0(3)=6000-V2; % <=0 c0(4)=V2-18000; % <=0 c0(5)=6000-V3; % <=0 c0(6)=V3-18000; % <=0 c0(7)=6000-V4; % <=0 c0(8)=V4-18000; % <=0 c0(9)=6000-V5; % <=0 c0(10)=V5-18000;% <=0 c0(11)=10000-V6; % <=0 No es necesario modelar una restricción de igualdad %Creamos la matriz "c" con "ceros" en caso cumpla las restricciones y "unos" en caso viole las restricciones. for i=1:length(c0) if c0(i)>0 c(i)=1; else c(i)=0; end end %Penalizamos las violaciones a las restricciones. penalicacion=100000; %Función Objetivo de la posición "x" f=of+penalicacion*sum(c); TÓPICOS EN SISTEMAS EXPERTOS E INTELIGENCIA ARTIFICIAL function f=Resultados(x) %Graficando reusltados de generación D=[600 1000 900 500 400 300]; %Actualizar la demanda for i=1:length(x) gen(i,1)=x(i); end for e=1:length(x) gen(e,2)=D(e)-x(e); end bar(gen,'stack') title('Resumen de Resultados - Generación') xlabel('Periodo "j"') ylabel('(MW)') legend('Planta de Vapor (MW)','Planta Hidroeléctrica (MW)') grid on %Graficando Volúmenes V(1)=10000; %Calculamos los volúmenes. for i=1:length(x) if gen(i,2)==0 V(i+1)=V(i)+4*(1000); %q1=0; else V(i+1)=V(i)+4*(1000-(260+10*gen(i,2))); %q1=260+10*Ph1; end end figure plot(V) title('Resumen de Resultados - Volumen(j)') xlabel('Periodo "j"') ylabel('Volúmen (acre-ft)') grid on