optim_control.html

> restart;

> with(CodeGeneration):

> with(linalg);

-------------------------------------------------------------------------------------------------------------

Exerc�cio de Otimiza��o Din�mica

Giovani Tonel ([email protected]) - Outubro de 2006

--------------------------------------------------------------------------------------------------------------

Equa��es do Balan�o de Massa e Energia

> eqf1:=Diff(h,t)=(Fe/Area-Fs);
eqf2:=Diff(Ca,t)=(1/tau)*(Cae-Ca)-k*Ca;
eqf3:=Diff(T,t)=(1/tau)*(Te-T)+(-Hr)*k*Ca*(1/(ro*Cp))-U*At*(T-Tw)*(1/(ro*Area*h*Cp));

> tau:=Area*h/Fe;
k:=ko*exp(-Ea/(R*T));

Equa��es do Controlador

> Fs:=(Cv/Area)*(1-u1)*(h)^0.5;
Tw:=u2*(Tmax-Tmin)+Tmin;

Valores das variav�is do estado est�cion�rio:

> estac:={h=hss,T=Tss, Ca=Cass, Tw=Twss, u1=u1ss, u2=u2ss};

> u2ss:=Twss;

> Fss=subs(estac,eval(Fs));

> eqf1:=Diff(h,t)=(Fe/Area-Fs);
eqf2:=Diff(Ca,t)=(1/tau)*(Cae-Ca)-k*Ca;
eqf3:=Diff(T,t)=(1/tau)*(Te-T)+(-Hr)*k*Ca*(1/(ro*Cp))-U*At*(T-Tw)*(1/(ro*Area*h*Cp));

> eq5:=Diff(Y,T)=F(Y,U);
Y=matrix([[H],[CA],[T]]);
U=matrix([[u1],[u2]]);

> f:=vector([rhs(eqf1),rhs(eqf2),rhs(eqf3)]);

Vari�veis de Estado (x):

> y:=vector([h,Ca,T]);

Variaveis Manipuladas u:

> u:=vector([u1,u2]);

Aplicando Jacobiano no vetor f em rela��o a y,u, temos:

> jacobian([f1(y,u),f2(y,u),f3(y,u)],[y,u]);

Aplicando Jacobiano no vetor f em rela��o as variaveis de estado, temos:

> jacfy:=jacobian(f,y);

Substituindo os valores das vari�veis no estado estacion�rio no Jacobiano de f, temos:

> A:=subs(estac,evalm(jacfy));

Aplicando Jacobiano no vetor f em rela��o as variaveis manipuladas, temos:

> jacfu:=jacobian(f,u);

Substituindo os valores das vari�veis no estado estacion�rio no Jacobiano de f em rela��o as manipuladas, temos:

> B:=subs(estac,evalm(jacfu));

Exportando para o MatLab

> Fortran(A,resultname="A");

     A(1,1) = -0.5D0 * dble(Cv) / dble(Area) * dble(1 - u1ss) * dble(hs
    #s ** (-0.5D0))
     A(2,1) = -1 / Area / hss ** 2 * Fe * (Cae - Cass)
     A(2,2) = -dble(1 / Area / hss * Fe) - ko * exp(-Ea / R / Tss)
     A(2,3) = -ko * Ea / R / Tss ** 2 * exp(-Ea / R / Tss) * dble(Cass)
     A(3,1) = -0.1D1 / dble(Area) / dble(hss ** 2) * dble(Fe) * (Te - T
    #ss) + U * At * (Tss - Twss * (Tmax - Tmin) - Tmin) / ro / dble(Are
    #a) / dble(hss ** 2) / Cp
     A(3,2) = -Hr * ko * exp(-Ea / R / Tss) / ro / Cp
     A(3,3) = -dble(1 / Area / hss * Fe) - Hr * ko * Ea / R / Tss ** 2
    #* exp(-Ea / R / Tss) * dble(Cass) / ro / Cp - U * At / ro / dble(A
    #rea) / dble(hss) / Cp