/home/ask/ask/ICT-Chemengg-resources-11062024/Sem6/Chemical-Process-Control/python/GP T-assisted-CSTR-PID-total_controller-FINAL.

py
Page 1 of 6 Sunday 02 February 2025 08:21:51 PM

1
2 import numpy as np
3 from scipy.integrate import solve_ivp
4 from scipy.optimize import fsolve
5 import matplotlib.pyplot as plt
6 # Reactor and physical constants
7 F = 100 # Flow rate [L/min]
8 Cai = 1 # Feed concentration [mol/L]
9 Ti = 273.15 + 80 # Inlet temperature [K]
10 V = 100 # Reactor volume [L]
11 # Reaction kinetics and thermodynamics
12 ko = 7.2e11 # Reaction rate constant pre-exponential factor [min^-1]
13 E = 55000 # Activation energy [J/mol]
14 R = 8.314 # Universal gas constant [J/(mol*K)]
15 delH = -6572.343578 * 4.184 # Reaction enthalpy [J/mol] (converted to J/mol)
16
17 # Heat transfer and physical properties for reactor and coolant
18 rho = 1 # Density [kg/L]
19 Cp = 1000 * 4.1842 # Heat capacity [J/(kg*K)]
20 UA = 7000 * 5.75 * 4.1842 # Overall heat transfer coefficient * area [J/(min*K)]
21 # Coolant properties
22 Fc = 15 # Coolant flow rate [L/min]
23 Vj = 15 # Jacket volume [L]
24 Tci = 273.15 + 18 # Coolant inlet temperature [K]
25 rhoc = rho # Coolant density [kg/L] (assumed same as reactor fluid)
26 Cpc = Cp # Coolant heat capacity [J/(kg*K)]
27 def steady_state_equations(vars):
28 Cass, Tss, Tcss = vars
29
30 # Reaction rate term
31 reaction_rate = ko * Cass * np.exp(-E / (R * Tss))
32
33 # Equation 1: Species balance in the reactor
34 eq1 = (F * (Cai - Cass) / V) - reaction_rate
35
36 # Equation 2: Energy balance in the reactor
37 eq2 = (F * (Ti - Tss) / V) - ((reaction_rate * delH) / (rho * Cp)) - (UA *
(Tss - Tcss) / (rho * V * Cp))

- 1 -
/home/ask/ask/ICT-Chemengg-resources-11062024/Sem6/Chemical-Process-Control/python/GP T-assisted-CSTR-PID-total_controller-FINAL.py
Page 2 of 6 Sunday 02 February 2025 08:21:51 PM

38
39 # Equation 3: Energy balance in the coolant jacket
40 eq3 = (Fc * (Tci - Tcss) / Vj) + (UA * (Tss - Tcss) / (rhoc * Cpc * Vj))
41
42 return [eq1, eq2, eq3]
43 # Initial guess for steady-state [Cass, Tss, Tcss]
44 initial_guess = (0.5, 350, 300)
45 Cass, Tss, Tcss = fsolve(steady_state_equations, initial_guess)
46 print("Steady-state solution:")
47 print(f" Cass = {Cass:.6f} mol/L")
48 print(f" Tss = {Tss:.4f} K")
49 print(f" Tcss = {Tcss:.4f} K")
50 # Perturbation values for simulation (setpoint changes)
51 Cai_shift = 0
52 Ti_shift = 0
53 Tci_shift = 2
54 T_setP = -3
55 # Controller parameters for reactor and coolant flows
56 controller_params = {
57 'kcF': 1, # Proportional gain for reactor flow
58 'tauiF': 1000, # Integral time for reactor flow
59 'taudF': 0, # Derivative time for reactor flow
60 'kvF': 6.6666, # Gain factor for reactor flow
61
62 'kcFc': -1, # Proportional gain for coolant flow
63 'tauiFc': 1000, # Integral time for coolant flow
64 'taudFc': 0, # Derivative time for coolant flow
65 'kvFc': 2.5, # Gain factor for coolant flow
66 }
67 # Global dictionary to store controller adjustments over time
68 adjustments_data = {
69 't': [],
70 'F': [],
71 'Fc': []
72 }
73 class Controller:
74
75 def __init__(self, params):

- 2 -
/home/ask/ask/ICT-Chemengg-resources-11062024/Sem6/Chemical-Process-Control/python/GP T-assisted-CSTR-PID-total_controller-FINAL.py
Page 3 of 6 Sunday 02 February 2025 08:21:51 PM

76 self.kcF = params['kcF']
77 self.tauiF = params['tauiF']
78 self.taudF = params['taudF']
79 self.kvF = params['kvF']
80
81 self.kcFc = params['kcFc']
82 self.tauiFc = params['tauiFc']
83 self.taudFc = params['taudFc']
84 self.kvFc = params['kvFc']
85
86 self.integral_error = 0.0
87 self.prev_error = T_setP # initialize previous error
88
89 def update(self, error, dt):
90 # Update integral term
91 self.integral_error += error * dt
92
93 # Calculate derivative term
94 derivative_error = (error - self.prev_error) / dt
95 self.prev_error = error
96
97 # Calculate raw adjustments (controller outputs) for reactor and coolant
flows
98 Gc_F = self.kcF * (error + (self.integral_error / self.tauiF) +
(self.taudF * derivative_error))
99 Gc_Fc = self.kcFc * (error + (self.integral_error / self.tauiFc) +
(self.taudFc * derivative_error))
100 # Apply saturation limits
101 # Gc_F which is differential P deviation: must be within [-3, 3]
102 Gc_F = max(min(Gc_F, 3), -3)
103 # Gc_Fc which is differential P deviation: must be within [-3, 3]
104 Gc_Fc = max(min(Gc_Fc, 3), -3)
105 F_adjustment = self.kvF * Gc_F
106 Fc_adjustment = self.kvFc * Gc_Fc
107 return F_adjustment, Fc_adjustment
108
109
110 def cstr_jacketed_odes(t, sys, controller, dt):

- 3 -
/home/ask/ask/ICT-Chemengg-resources-11062024/Sem6/Chemical-Process-Control/python/GP T-assisted-CSTR-PID-total_controller-FINAL.py
Page 4 of 6 Sunday 02 February 2025 08:21:51 PM

111 Ca_dev, T_dev, Tc_dev = sys

112
113 # Compute control error (using reactor temperature deviation)
114 error = T_setP - T_dev
115
116 # Get control adjustments for flows
117 F_adjust, Fc_adjust = controller.update(error, dt)
118
119 # Record the controller adjustments along with time
120 adjustments_data['t'].append(t)
121 adjustments_data['F'].append(F_adjust)
122 adjustments_data['Fc'].append(Fc_adjust)
123
124 # Reaction rate based on steady state values (assuming deviations are small)
125 #reaction_rate_ss = ko * Cass * np.exp(-E / (R * Tss))
126
127 # Mass balance: deviation in concentration
128 dCa_dt = ((F*Cai_shift)/V) + (((Cai-Cass)*F_adjust)/V) - (Ca_dev*((F/V) +
(ko*np.exp(-E/(R*Tss))))) -
(T_dev*((Cass*E*ko*np.exp(-E/(R*Tss)))/(R*(Tss**2))))
129
130 # Energy balance for reactor (deviation dynamics)
131 dT_dt = ((F_adjust*(Ti-Tss))/V) + ((Ti_shift*F)/V) + ((Tc_dev*UA)/(rho*V*Cp))
- (Ca_dev*((delH*ko*np.exp(-E/(R*Tss)))/(rho*Cp))) -
(T_dev*((UA/(rho*V*Cp))+(F/V)+((delH*E*ko*Cass*np.exp(-E/(R*Tss)))/(rho*Cp*R*(Ts
s**2)))))
132
133 # Energy balance for coolant (deviation dynamics)
134 dTc_dt = ((Fc_adjust*(Tci-Tcss))/Vj) + ((Tci_shift*Fc)/Vj) +
((T_dev*UA)/(rhoc*Cpc*Vj)) - (Tc_dev*((Fc/Vj)+(UA/(rhoc*Cpc*Vj))))
135
136 return [dCa_dt, dT_dt, dTc_dt]
137
138 t_final = 18
139 num_points = 500
140
141 def run_simulation(t_final, num_points):
142 t_eval = np.linspace(0, t_final, num_points)

- 4 -
/home/ask/ask/ICT-Chemengg-resources-11062024/Sem6/Chemical-Process-Control/python/GP T-assisted-CSTR-PID-total_controller-FINAL.py
Page 5 of 6 Sunday 02 February 2025 08:21:51 PM

143 dt = t_eval[1] - t_eval[0]

144
145 # Reset the global adjustments_data dictionary in case of multiple runs
146 adjustments_data['t'].clear()
147 adjustments_data['F'].clear()
148 adjustments_data['Fc'].clear()
149
150 # Initial deviations from steady state (zero deviations initially)
151 Ca0_dev = 0.0
152 T0_dev = 0.0
153 Tc0_dev = 0.0
154 initial_state = [Ca0_dev, T0_dev, Tc0_dev]
155
156 # Instantiate controller
157 controller = Controller(controller_params)
158
159 # ODE solver using RK45 method
160 sol = solve_ivp(fun=lambda t, y: cstr_jacketed_odes(t, y, controller, dt),
161 t_span=(t_eval[0], t_eval[-1]),
162 y0=initial_state,
163 t_eval=t_eval,
164 method='RK45')
165
166 # Recover deviations and add back the steady-state values
167 Ca_total = Cass + sol.y[0] # Total concentration signal [mol/L]
168 T_total = Tss + sol.y[1] # Total reactor temperature [K]
169 Tc_total = Tcss + sol.y[2] # Total coolant temperature [K]
170
171 return sol.t, Ca_total, T_total, Tc_total
172
173
174 def plot_results(t, Ca, T, Tc):
175 plt.figure(figsize=(10, 8))
176 # Plot reactor concentration
177 plt.subplot(3, 1, 1)
178 plt.plot(t, Ca, label=f"Reactor Concentration (Cass = {Cass:.6f} mol/L)")
179 plt.xlabel("Time [min]")
180 plt.ylabel("Concentration [mol/L]")

- 5 -
/home/ask/ask/ICT-Chemengg-resources-11062024/Sem6/Chemical-Process-Control/python/GP T-assisted-CSTR-PID-total_controller-FINAL.py
Page 6 of 6 Sunday 02 February 2025 08:21:51 PM

181 plt.legend()
182 plt.grid(True)
183
184 # Plot reactor and coolant temperatures
185 plt.subplot(3, 1, 2)
186 plt.plot(t, T, label=f"Reactor Temperature (Tss = {Tss:.3f} K)")
187 plt.plot(t, Tc, label=f"Coolant Temperature (Tcss = {Tcss:.3f} K)")
188 plt.xlabel("Time [min]")
189 plt.ylabel("Temperature [K]")
190 plt.legend()
191 plt.grid(True)
192
193 # Plot controller adjustments
194 plt.subplot(3, 1, 3)
195 plt.plot(adjustments_data['t'], adjustments_data['F'], label="F_adjustment")
196 plt.plot(adjustments_data['t'], adjustments_data['Fc'], label="Fc_adjustment")
197 plt.xlabel("Time [min]")
198 plt.ylabel("Adjustment Value")
199 plt.legend()
200 plt.grid(True)
201
202 plt.tight_layout()
203 plt.show()
204
205 def main():
206 # Run dynamic simulation; adjust t_final as needed (e.g., 30 minutes)
207 t, Ca, T, Tc = run_simulation(t_final, num_points)
208 # Plot the simulation results including controller adjustments
209 plot_results(t, Ca, T, Tc)
210
211 if __name__ == '__main__':
212 main()
213

- 6 -