Transcript Aplicações em MATLAB e SIMULINK

MATLAB EM SIMULAÇÃO E
CONTROLE DE PROCESSOS
PROFa. OFÉLIA DE Q.F. ARAÚJO
[email protected]
www.eq.ufrj.br/links/h2cin/eqe487
Desafio
Solução
Fenômeno
Físico/
Químico
Programação
Métodos
Numéricos
Desafio
Dinâmica
Linear e
Não-Linear
Dinâmica
Linear
Provas e
Trabalhos
em GRUPO
Avaliação
Exclusivamente com
PROVAS
PASCAL
Métodos
Numéricos
MATLAB
Built in
Functions
(ODE Suite)
Desafio
Controle de
Processos:
Solução e
SIMULINK
Toolbox de
Controle
Modelagem
Fenômeno
de Processos:
Físico/
MATLAB e
Químico
SIMULINK
Intro. Eng.
Quim.:
ProgramaEnsino
ção
Programar
c/ MATLAB
Rotinas
Métodos
MATLAB:
ODE Suite,
Numéricos
etc
Caso I
Ensino de
Programação
Introdução
a Cálculos
em Eng.
Quím.
Experimentos
COMPUTA
-CIONAIS
em EQ
Introdução à Engenharia Química:
rápido aprendizado na solução
de problemas introdutórios
Caso I
1
• Ensino de Programação com MATLAB por EQ.
• Substitui PASCAL
2
• Ensino de Programação com Exemplos de EQ.
• Resolução de Balanços de Massa – Equações Lineares
3
• Conteúdo de Engenharia Química no início do Curso
• Redução da Evasão
Caso I
Modelagem & Dinâmica de P
Caso I
Modelos simples - o tanque de
pode-se escrever o balanço de massa do sistema
dmt 
  FE  F 
dt
Ainda,
dmt 
dht 
 A
dt
dt
e, portanto,
dht  1
 FE  F 
dt
A
Caso I
Modelagem & Dinâmica de Processos
Tanque vazio no início do
experimento,
Modelos simples - o tanque de nível
Freqüentemente, considera-se a vazão de saída do tanque
proporcional
à
altura
da
coluna
de
líquido
é inversamente
quando passa a ser alimentado
Modelagem
com vazão FE & Dinâmica de
proporcional a uma resistência ao escoamento (R):
h
F 
R
Logo,
FE é constante
Modelos simples - o tanque d
(5)
Este modelo simples de um tanque de ní
energia, possui uma solução analítica:
dht  1 
h
  FE  
dt
A
R
Solução(6)
Analítica
1
t



ht   RFE 1  e RA 


Para simular este modelo, basta escol
constantes R, A e FE, das condições iniciais h0 e
A simulação da solução analítica do mo
nível é mostrada a seguir.
Modelagem & Dinâmica de Process
Modelos simples - o tanque de nível
Caso I
Este modelo simples de um tanque de nível, sem bala
% Definição das constantes do modelo
energia, possui uma solução analítica:
R = 1;
% h/m2
t



RA
ht   RFE 1  e 
A = 2;
% m2


Fe = 10;
% m3/h
% Tempo de simulação
Para simular este modelo, basta escolher os valore
t = 0.0:0.01:10.0;
% constantes
h
R, A e FE, das condições iniciais h0 e t0.
A simulação da solução analítica do modelo do tanq
% Simulação da altura de líquido
nível é mostrada a seguir.
h = R*Fe*(1 - exp(-t/(R*A)));
% m
% Visualização da simulação
plot(t,h)
title('Simulação do tanque de nível')
xlabel('Tempo (h)')
ylabel('Altura (m)')
Caso I
Modelagem & Dinâmica de Processos
Modelos simples - o tanque de nível
Este modelo simples de um tanque de nível, sem balanço de
energia, possui uma solução analítica:
t



ht   RFE 1  e RA 


(7)
Para simular este modelo, basta escolher os valores das
constantes R, A e FE, das condições iniciais h0 e t0.
A simulação da solução analítica do modelo do tanque de
nível é mostrada a seguir.
1
Caso II
1
2
3
• Apresentação de Sistemas Dinâmicos
• Aplicação de Leis de Conservação (Massa, Energia e Momento)
• Introdução das rotinas da ODE Suite e SIMULINK
• Resolução de Balanços – EDO’s Não-Lineares Acopladas
• Simulação de Sistemas Não-Lineares com Mínimo Esforço de
Programação
• Experimentos com dinâmica de processos (instabilidades, ciclos, etc)
Caso II
Fenôm.
Físico/
Químicos
Leis de
Conservação
Experimentos de
DINÂMICA
em EQ
Modelagem e Dinâmica de Processos: como
rapidamente um aluno pode simular um
processo, avaliando comportamentos
dinâmicos complexos.
Caso II
Caso IIa
2 Tanques em série, volumes V1 e V2, vazão volumétrica F e concentração C
dC
1  FC  FC  (1   )FC
V
1 dt
2
0
1
F
dC
2  (1   ) F (C (t )  C (t ))
V
2 dt
1
2
F
C (t)
0
V
1
C
1
V
2
C (t)
1
C
2
F
C (t)
2
Caso IIa
Caso IIa
demo
Caso IIa
Cˆ ( s)
K
G( s)  2

Cˆ ( s) ( 1s  1)( 2 s  1)  K
0
Caso IIa
xsi=[0.2 0.6 0.9 1. 2]
K=2;tau=10;
for i=1:length(xsi)
G(i)=tf([K],[tau^2 2*xsi(i)*tau 1]);
end
step(G)
Caso IIb
dV
 (FE  F )
dt
dCR
E
 (FECR ,E  CR F ) / V  k0 exp(
)CR
dt
RT
V
dTR
E


 (FETR,E  TFE )  Hr k0 exp(
)CR  UA(T  Th ) / Cp
dt
RT


Caso IIb
Parâmetros
function [sys,x0] = reator(t,x,u,flag,U,A,DeltaH,ro,Cp,E,R,k0)
%
% Simula um reator CSTR (mistura perfeita) no qual se conduz uma
% reação exotérmica (A->B), resfriado por serpentina
%
% U
= 150
BTU/(h.ft2.R), coeficiente de troca térmica
% A
= 250
ft2,
área de troca térmica
% DeltaH = -30000 BTU/lbm,
calor de reação
% ro
= 50
lb/ft3,
densidade
% Cp
= 0.75
BTU/(lbm.R),
calor específico
Blá-blá-blá
% E
= 30000
BTU/lbm,
energia de ativação
% R
= 1.99
BTU/(lbm.R),
constante dos gases
% k0
= 7.08e10 1/h,
termo pré-exponencial da constante de reação
%
switch flag
case 0 % Dimensiona o sistema e inicializa os estados
Condições
iniciais
% sys=[estados,0,saídas,entradas,0,0]
sys = [3,0,3,5,0,0];
% Condições iniciais
cr = 0.1315;
%lbm/ft3, concentração inicial no reator
T
= 500;
%R,
temperatura do reator
V
= 200;
%ft3,
volume do reator
x0 = [cr T V];
Caso IIb
dV
 (FE  F )
dt
end
Entradas
case 1 % Calcula as derivadas
% Atualiza entradas
Cre = u(1);
%lbm/ft3, concentração da alimentação
Fe = u(2);
%ft3/hr, vazão de alimentação
F
= u(3);
%vazão de retirada
Tc = u(4);
%R, temperatura do fluido de refrigeração
Te = u(5);
%R, temperatura da alimentação
% Cálculo das derivadas
O
Cr = x(1);
MODELO!!
T = x(2);
V = x(3);
dCR
E
k = k0*exp(-E/(R*T));
 (FECR ,E  CR F ) / V  k0 exp(
)CR
dt
RT
dCr = (Fe*Cre-F*Cr)/V) - k*Cr;
dV = Fe-F;
dT = (Fe*Cp*ro*(Te-T) + DeltaH*k*Cr*V - U*A*(T-Tc)) /(V*ro*Cp);
sys = [dCr; dT; dV];
dT
E


V R  (FETR,E  TFE )  Hr k0 exp(
)CR  UA(T  Th ) / Cp
case 3 % Calcula as saídas
dt
RT


sys = [x(1) x(2) x(3)];
otherwise
sys = [];
Reator.m
demo
Caso IIc
U s 
G1 s
G2 s
2
G1 s  
5s  1


Y s 
1
G2 s  
s 1
demo
Caso IIc
G1=tf(2,[5 1]);
G2=tf(-1,[1 1]);
step(G1+G2)
2
G1 s  
5s  1
1
G2 s  
s 1
Caso IId
Caso IId
“Experimento” para geracão de dados experimentais:
Temperatura do Tanque 1
Temperatura do Tanquer 2
%Roda Laboratorio de EQ
sim('ldeq');
dados=[T1 T2(:,2)];
plot(T1(:,1),T1(:,2),'*b'),
xlabel('Tempo'), ylabel('T1')
save dados dados
demo
Caso IId
demo
Caso IId
%Obtenção de Funções de Transferência a partir de respostas
%temporais a perturbação degrau.
%Carrega séries temporais
load dados
tempo=dados(:,1);
y=dados(:,2);
degrau=2;
%Valor inicial para parâmetros
Tau=10;
Ganho=10;
Tmorto=10;
X0=[Tau Ganho Tmorto];
opcao=optimset('Display','iter');
Xotimo=fminsearch('fobj',X0,opcao,degrau,y,tempo)
Caso IId
function [f] = fobj(X,degrau,y,tempo)
Tau
Tempo
ycalc
= abs(X(1)); Ganho = abs(X(2)); TM
= max(0,(tempo-TM));
= degrau*Ganho*(1-exp(-tempo/Tau));
= abs(X(3));
%Somatório do quadrado dos erros
f=(y-ycalc)'*(y-ycalc);
%DESENHA PARA ACOMPANHAR A OTIMIZAÇÃO
plot(tempo,y,'*',tempo,ycalc,'m')
%ESCREVE NO GRÁFICO OS VALORES DOS PARÂMETROS
title(num2str(X))
%A rotina a seguir força o MATLAb a desenhar o gráfico imediatamente
drawnow
Caso IId
Caso III
Modelo
Dinâmico
Lei de
Controle
Experimentos de
CONTROLE
em EQ
Dinâmica e Controle de Processos: simular
estratégias de controle e sintonizar
controladores com MATLAB, SIMULINK e
Toolboxes de Controle e Otimização.
Caso III
1
2
3
• Análise de sistemas dinâmicos
• Proposição de algoritmos de controle
• Introdução ao Toolbox de Controle
• Análise da respostas dinâmicas de malhas de controle
• MATLAB/SIMULINK
• Experimentos com sintonia de controladores , otimização de indicadores
de desempenho de malhas de controle
Caso III
Caso IIIa
Expansão em Frações Parciais ....

R3
1
s2  s  1
R1
R2
R4
X (s )   3





s  s  9s2  26s  24  s  P1 s  P2 s  P3 s  P4
• Expansão em
Frações
Parciais
1
2
• Inversão para
o Domínio do
Tempo
• Simular para
horizonte de
resposta
3
Caso IIIa
Expansão em Frações Parciais ....
R3
s2  s  1
R1
R2
R4
X (s )  4




s  9s3  26s2  24s s  P1 s  P2 s  P3 s  P4
% Coeficientes de P(s) e Q(s) em ordem decrescente de
potências de s
Den = [1 9 26 24 0];
Num = [1 1 1];
[R,P,K] = residue(Num,Den)
Caso IIIa
Expansão em Frações Parciais .... E inversão para o domínio do Tempo
R =
-1.6250
2.3333
-0.7500
0.0417
P =
-4.0000
-3.0000
-2.0000
0
K =
[]
X (s ) 
R3
R1
R2
R4



s  P1 s  P2 s  P3 s  P4
x(t )  1,650e4t  2,3333e3t  0,75e2t  0,0417
Caso IIIa
% Simulação
Den = [1 9 26 24 0];
Num = [1 1 1];
[R,P,K] = residue(Num,Den)
Tempo=0:0.1:10;
X=zeros(size(Tempo));
for i=1:length(R)
X=X+R(i)*exp(P(i)*Tempo);
end
plot(Tempo,X)
Caso IIIa
% Estabilidade
Den = [1 9 26 24 0];
Num = [1 1 1];
roots(Den)
ans =
0
-4.0000
-3.0000
-2.0000
Caso IIIa

1
s2  s  1
X (s )   3

s  s  9s2  26s  24 
• Diagrama de
Blocos
1
2
• SIMULINK
• Simular para
horizonte de
resposta
3
Caso IIIa
1
s
s2  s  1
s3  9s2  26s  24
U(S)
Step
s2 +s+1
s3 +9s2 +26s+24
Transfer Fcn
demo
X(S)
Scope
Caso IIIa
Step
s2 +s+1
s3 +9s2 +26s+24
Transfer Fcn
Scope
Caso IIIb
Estabilidade de Malha de Controle
demo
kc=[3 4 5 6 7 8];
simb=['c', 'm', 'g', 'b', 'r', 'k'];
Texto=[];
kp=8;
den=conv(conv([1 2],[1 2]),[1,2]);
num=kp;
Gp=tf(num,den);
Gv=tf(1,1);
Gm=tf(1,1);
t=linspace(1,10,100);
hold on
for i=1:length(kc)
Gc=tf(kc(i),1);
G=Gc*Gv*Gp/(1+Gc*Gv*Gp*Gm);
[y,t]=step(G,t);
Texto=[Texto; ['K_C = ' num2str(kc(i))]];
plot(t,y,simb(i),'LineWidth',2)
end
legend(Texto)
xlabel('Tempo')
ylabel('Variável Controlada')
Caso IIIb
Estabilidade de Malha de Controle
demo
G( s ) 
Equação
Característica
Gc ( s)Gv ( s)Gp ( s)
1  Gc ( s)Gv ( s)Gp ( s)Gm ( s)
tau = 0.5;
cor=['m';'g';'k';'c';'y';'b';'r'];
cor=[cor; cor; cor];
xsi = 0 :0.2 :2;
Raizes=[];
for i=1:length(xsi)
% Escreve-se a equação característica como: tau^2 + 2*xsi*tau +1
Raizes = [Raizes roots([tau^2, 2*xsi(i)*tau, 1])];
h=plot(real(Raizes(:,i)),imag(Raizes(:,i)),[cor(i,:) 's']);hold on
set(h,'Markersize',10,'Markerfacecolor',cor(i,:),'Markeredgecolor','k');
legenda{i}=['xsi = ' num2str(xsi(i))];
end
h=legend(legenda);
h=plot(real(Raizes(1,:)),imag(Raizes(1,:)));
set(h,'Linestyle','-','LineWidth',2);
h=plot(real(Raizes(2,:)),imag(Raizes(2,:)));
set(h,'Linestyle','-','LineWidth',2);
xlabel('Real');
ylabel('Imaginário');
grid on
G( s ) 
Gc ( s)Gv ( s)Gp ( s)
1  Gc ( s)Gv ( s)Gp ( s)Gm ( s)
xsi=[0.3 0.6 0.8 1 1.2];
tau = 2;
for i=1:length(xsi)
den=[tau^2 2*xsi(i)*tau 1];
G(i)=tf([1],den);
end
bode(G(1),G(2),G(3),G(4),G(5))
nyquist(G(1),G(2),G(3),G(4),G(5))
Demonstrações
demo
Demonstrações
demo
Sistemas Multicapacitivos
% Multicapacitivos
g1=tf(1,[10 1])
g2=g1*g1
g3=g1*g1*g1
g4=g1*g1*g1*g1
step(g1,g2,g3,g4)
legend({'g1'; 'g2'; 'g3'; 'g4'})
Sintonia de PID
Neste exemplo, aplica-se procedimento de
sintonia baseado em otimização do critério de
desempenho ISE. Utiliza-se modelo da malha
feedback em ambiente SIMULINK para
processo com resposta inversa descrito por:
(1  3 s)
G( s) 
(2s  1)(5 s  1)
O denominador é reescrito para o formato polinomial, obtendo-se os
coeficientes com a rotina de convolução (conv) do MATLAB:
>> conv([2 1], [5 1])
ans =
10
7
1
demo
otimiza.m
warning off
global Kc TauI TauD
% Abre o modelo do SIMULINK
pid
% Atribui valores iniciais aos parâmetros de sintonia
Kc
= 1;
TauI = 1;
TauD = 1;
teta = [Kc; TauI; TauD];
% Especifica a janela de otimização (tempo de simulação)
tfinal=100;
% Configura parâmetros do otimizador e chama otimizador
options=optimset('MaxFunEvals',200,'Display','iter');
[tetaotim] = fminsearch('fob', teta, options, tfinal)
fob.m
function f
= fob(teta, tfinal)
global Kc TauI TauD
% Recebe o novo conjunto de parametros
Kc
= abs(teta(1));
TauI = abs(teta(2));
TauD = abs(teta(3));
% Simula com os novos parametros
[t,x,y] = sim('pid',tfinal);
% Processa os resultados do após a simulação:
% Integração empregando área de retângulo
SE
= y.*y;
ISE = sum(SE(1:end-1).*(t(2:end)-t(1:end-1)));
f=ISE;
demo
Parâmetros
Kc
I
D
Valor Inicial
-1,0
1,0
1,0
Valor Otimizado
-11,2
8,1
2,5
Demonstrações – Controle Cascata
demo