kill(all)$

Temperature control of a reactor equipment , 'version 1.0'
peter.vlasschaert@gmail.com,21/12/2019

load packages
---------------------

load(lrats)$

plant
====

Ti  = inlet temperature
Tis = assumed inlet temperature
ρ    = density ( = kg/m^3)
 V  = volume  ( = m^3)
C  = heat capacity ( = J/(kg*K))
w  = kg/s
qs  = heat stream  ( = J/s)

control
=====

TR(t) = reference temperature for T by negative feedback

We want to solve differential equations with 'Laplace Transformation =LT'
' independent variable must be declared explicit '   T = > T(t) , q => q(t)

1e)  odes
------------

plant unit
------------

p1:ρ*V*C*diff(T(t),t) = w*C*(Ti-T(t))+q(t);

control unit 
---------------

p1a:q(t) =w*C*(TR(t)-Tis)+Kc*(TR(t)-T(t));

2e) conversion  ( to impliment time constant)
-------------------------------------------------------------

plant 
====

p2:p1+w*C*T(t),expand;

p2a:p2/( C*w),expand;

use : time constant : τ1 = V*ρ/w

p2b:lratsubst([ V*ρ/w =τ1],lhs(p2a)) =rhs(p2a);

3e) LT
----------

plant
====

p3:laplace(p2b,t,s);

suppose : T(0) =0

p3a:factor(subst([T(0)=0],lhs(p3)) )=rhs(p3);

control
=====

p3b1:qs(s) =laplace(part(rhs(p1a),1),t,s);

p3b2:laplace(lhs(p1a),t,s) = lhs(p3b1) +laplace(part(rhs(p1a),2),t,s);

combine : plant & control, ' p3a,p3b2'
============================

p4:lratsubst([p3b2],p3a);

p4a:lhs(p4)= expand(rhs(p4));

p4b:p4a-laplace(T(t),t,s);

p4c:expand(p4b-part(rhs(p4b),2)-part(rhs(p4b),4));

p4d:part(lhs(p4c),1)+factor(part(lhs(p4c),2)+part(lhs(p4c),3))=rhs(p4c);

change ; to combination 'control & plant'
===============================

a)   qs(s) = 0 ,put disturbance to zero.
---------------------------------------------------

p5:subst([qs(s) = 0],p4d);

b)  (Tis = Tis+ΔTi,TR(s) = TR(s) ) ,@t = 0
-----------------------------------------------------------
result : 'p5a' substitution

p5a:subst([ Ti = ΔTi,laplace(TR(t),t,s) = 0],p5);

p5b:solve(p5a,laplace(T(t),t,s))[1];

c) additiv principal  ( = linear systems)
:=============================
input : {u1,u2}  =>  output :{y1,y2}
result:
u1+u2  = > y1 +y2

d) homogenity princpal ( = linear systems)
==================
c = constant
input:{u} => output:{y}
result :
c*u  => c*y

e)  combination : c) & d)  => superposition principal
======================================

d) find from : p5b: (Inverse Laplace Transform:ilt, 'Maxima')
--------------------------------------------------------------------------------
@ t = 0
======

p5c1:ilt(p5b,s,t);

find limit : t-> infinity
----------------------------

p5c1a:limit(p5c1,t,inf);   

assume(C*w*(C*w+Kc)*τ1> 0);

otherwise : Maxima ask question for this  .

p5c1b:lhs(p5c1a) =ev(rhs(p5c1a),nouns);

@t ≠ 0  (use :p4d)
=====

rem  s    = > d/dt
         1/s = >  integrate

p5d1:s*laplace(T(t),t,s)*τ1+(laplace(T(t),t,s)*(C*w+Kc))/(C*w)=(Kc*laplace(TR(t),t,s))/(C*w)+qs(s)/(C*w)+Ti/s;

suppose:T(s)  =laplace(T(t),t,s)
                TR(s) =laplace(TR(t),t,s)

p5d2:[laplace(T(t),t,s)=T(s),laplace(TR(t),t,s)=TR(s)];

p5d3:lratsubst(p5d2,p5d1);

p5d4:expand(p5d3);

time domain :manuel
================

p6:τ1*diff(T(t),t) +(1+Kc/(C*w))*T(t) = Ti+qs(t)/(C*w) + (Kc*TR(t))/(C*w);

Transfer Function
==============

1e) suppose :qs(s) = 0  into  'p5d4'

2e)construct : TF(s) = Transfer Function = T(s)/TR(s) 

p7:subst([qs(s) = 0],p5d4);

p7a:solve(p7,T(s))[1];

We can construct  TF(s) when Ti = 0
--------------------------------------------------

p7b:subst([Ti = 0 ],p7a);

p7c:TF(s) =rhs(p7b)/TR(s);

time domain => ilt

p7d:ilt(p7c,s,t);

p7e: TF(t) =rhs(p7d);

when manuel : ode from 'p6'
----------------------------------------

p8:p6;

1e) qs(t) = 0

p8a:subst([qs(t)=0],p8);

superposition principal
-------------------------------

2e)    T (t)   = > Tstar(t)+ TRR
         TRR   => TRR
          Ti   =>  ΔTi + TRR
======================
TR(t) = constant here = TRR
         

find  :  Tstar(t),use 'p8a'

p9:[T(t) = Tstar(t)+TRR,Ti=ΔTi+TRR,TR(t)=TRR];
p9a:subst(p9,p8a);

p9b:ev(p9a,diff);

p9c:expand(lhs(p9b) - rhs(p9b))  =0;

p9d: p9c+ΔTi;

solve :  LT => ilt

1e)  LT
----------

p10:laplace(p9d,t,s);

a) IC : Initial Condition ,Tstar(0)=0
--------------------------------------------

p10a: subst([Tstar(0)=0],p10);

b) solve for 'laplace(Tstar(t),t,s) from 'p10a'
-----------------------------------------------------------

p10b:solve(p10a,'laplace(Tstar(t),t,s))[1];

2e)  ilt
----------

p10d:Tstar(t) = ilt(rhs(p10b),s,t);

final value theorema
================

p11:limit(p10d,t,inf);

how to do this  laplace domain:

p11a:limit(t);

definition :

p11a:limit(Tstar(t),t,inf) = limit (s*Tstar(s),s,0);

use :rhs  from p10b

p11b:lhs(p11a) = limit(s*rhs(p10b),s,0);

initial value theorema
================

p11c:limit(p10d,t,0);

final result
========

p12: T(t) =rhs(p10d)+TRR;

 ===========================================================