python-3.xgekko
How to add time delay in state space model
I recently try to control my TCLab Arduino through a state space model. I refer to State Space Models and transform a first order linear system (without time delay) into state space form, the control effect is pretty good. Now I want to use first order plus dead time model to control the lab, but I don't know how to transform FOPDT model into state space form. How to add time delay in state space model?
Here is the code:
import tclab
import numpy as np
import time
import matplotlib.pyplot as plt
from gekko import GEKKO
import json
# Connect to Arduino
# a = tclab.TCLab()
a = tclab.TCLabModel()
# Run time in minutes
run_time = 60.0
# Number of cycles
loops = int(60.0*run_time)
tm = np.zeros(loops)
# Temperature (K)
T1 = np.ones(loops) * a.T1 # temperature (degC)
Tsp1 = np.ones(loops) * 40.0 # set point (degC)
# heater values
Q1s = np.ones(loops) * 0.0
Q1_ss = 0
#########################################################
# Initialize Model
#########################################################
tau = 160.0
kp = 0.6
Am = np.zeros((1,1))
Bm = np.zeros((1,1))
Cm = np.zeros((1,1))
Am[0, 0] = - 1/tau
Bm[0, 0] = kp/tau
Cm[0, 0] = 1
# state space simulation
m = GEKKO(remote=False)
x,y,u = m.state_space(Am,Bm,Cm,D=None)
mv = u[0]
cv = y[0]
mv.VALUE = Q1_ss
mv.STATUS = 1 # use to control temperature
mv.FSTATUS = 0 # no feedback measurement
mv.LOWER = 0.0
mv.UPPER = 100.0
mv.DMAX = 10.0
mv.COST = 0.0
mv.DCOST = 0.1
cv.VALUE = a.T1
cv.STATUS = 1 # minimize error with setpoint range
cv.FSTATUS = 1 # receive measurement
cv.TR_INIT = 2 # reference trajectory
cv.TAU = 60 # time constant for response
m.time = np.linspace(0, 160, 81)
m.options.IMODE = 6 # MPC
m.options.CV_TYPE = 1 # Objective type
m.options.NODES = 2 # Collocation nodes
m.options.MAX_TIME = 10
##################################################################
# Create plot
plt.figure()
plt.ion()
plt.show()
filter_tc1 = []
def movefilter(predata, new, n):
if len(predata) < n:
predata.append(new)
else:
predata.pop(0)
predata.append(new)
return np.average(predata)
# Main Loop
start_time = time.time()
prev_time = start_time
last_Q1 = Q1_ss
try:
for i in range(1,loops):
# Sleep time
sleep_max = 2.0
print(time.time() - prev_time)
sleep = sleep_max - (time.time() - prev_time)
if sleep>=0.01:
time.sleep(sleep)
else:
time.sleep(0.01)
# Record time and change in time
t = time.time()
dt = t - prev_time
prev_time = t
tm[i] = t - start_time
# Read temperatures in Kelvin
curr_T1 = a.T1
last_T1 = curr_T1
avg_T1 = movefilter(filter_tc1, last_T1, 3)
T1[i] = curr_T1
###############################
### MPC CONTROLLER ###
###############################
cv.MEAS = avg_T1
# input setpoint with deadband +/- DT
DT = 0.5
cv.SPHI = Tsp1[i] + DT
cv.SPLO = Tsp1[i] - DT
# solve MPC
m.solve(disp=False)
# test for successful solution
if (m.options.APPSTATUS==1):
# retrieve the first Q value
Q1s[i] = mv.NEWVAL
print('Q1.NEWVAL', mv.NEWVAL)
with open(m.path+'//results.json') as f:
results = json.load(f)
else:
# not successful, set heater to zero
Q1s[i] = last_Q1
print('last_Q1', mv.NEWVAL)
print(m.path)
# Write output (0-100)
a.Q1(Q1s[i])
# Plot
plt.clf()
ax=plt.subplot(2,1,1)
ax.grid()
plt.plot(tm[0:i],T1[0:i],'ro',MarkerSize=3,label=r'$T_1$')
plt.plot(tm[0:i],Tsp1[0:i],'b-',MarkerSize=3,label=r'$T_1 Setpoint$')
cv_name = cv.NAME + '.bcv'
print(cv_name)
plt.plot(tm[i]+m.time,results[cv_name],'k-.',\
label=r'$T_1$ predicted',linewidth=3)
# plt.plot(tm[i]+m.time,results['tc1.tr_hi'],'k--',\
# label=r'$T_1$ trajectory')
# plt.plot(tm[i]+m.time,results['tc1.tr_lo'],'k--')
plt.ylabel('Temperature (degC)')
plt.legend(loc='best')
ax=plt.subplot(2,1,2)
ax.grid()
plt.plot(tm[0:i],Q1s[0:i],'r-',LineWidth=3,label=r'$Q_1$')
plt.plot(tm[i]+m.time,mv.value,'k-.',\
label=r'$Q_1$ plan',linewidth=3)
# plt.plot(tm[0:i],Q2s[0:i],'b:',LineWidth=3,label=r'$Q_2$')
plt.ylabel('Heaters')
plt.xlabel('Time (sec)')
plt.legend(loc='best')
plt.draw()
plt.pause(0.05)
# Turn off heaters
a.Q1(0)
a.Q2(0)
print('Shutting down')
a.close()
# Allow user to end loop with Ctrl-C
except KeyboardInterrupt:
# Disconnect from Arduino
a.Q1(0)
a.Q2(0)
print('Shutting down')
a.close()
# Make sure serial connection still closes when there's an error
except:
# Disconnect from Arduino
a.Q1(0)
a.Q2(0)
print('Error: Shutting down')
a.close()
Solution
Use the delay function in Gekko to add time delay. Here is an example. To add it to a state space model you can either delay the input mv or output cv. Here is the state space model with the delay on the output.
# state space simulation
m = GEKKO(remote=False)
x,y,u = m.state_space(Am,Bm,Cm,D=None)
mv = u[0]
cv_in = y[0]
cv = m.CV()
m.delay(cv_in,cv,4) # delay of 4 steps (8 sec)
import tclab
import numpy as np
import time
import matplotlib.pyplot as plt
from gekko import GEKKO
import json
# Connect to Arduino
# a = tclab.TCLab()
a = tclab.TCLabModel()
# Run time in minutes
run_time = 60.0
# Number of cycles
loops = int(60.0*run_time)
tm = np.zeros(loops)
# Temperature (K)
T1 = np.ones(loops) * a.T1 # temperature (degC)
Tsp1 = np.ones(loops) * 40.0 # set point (degC)
# heater values
Q1s = np.ones(loops) * 0.0
Q1_ss = 0
#########################################################
# Initialize Model
#########################################################
tau = 160.0
kp = 0.6
Am = np.zeros((1,1))
Bm = np.zeros((1,1))
Cm = np.zeros((1,1))
Am[0, 0] = - 1/tau
Bm[0, 0] = kp/tau
Cm[0, 0] = 1
# state space simulation
m = GEKKO(remote=False)
x,y,u = m.state_space(Am,Bm,Cm,D=None)
mv = u[0]
cv_in = y[0]
cv = m.CV()
m.delay(cv_in,cv,4) # delay of 4 steps (8 sec)
mv.VALUE = Q1_ss
mv.STATUS = 1 # use to control temperature
mv.FSTATUS = 0 # no feedback measurement
mv.LOWER = 0.0
mv.UPPER = 100.0
mv.DMAX = 10.0
mv.COST = 0.0
mv.DCOST = 0.1
cv.VALUE = a.T1
cv.STATUS = 1 # minimize error with setpoint range
cv.FSTATUS = 1 # receive measurement
cv.TR_INIT = 2 # reference trajectory
cv.TAU = 60 # time constant for response
m.time = np.linspace(0, 160, 81)
m.options.IMODE = 6 # MPC
m.options.CV_TYPE = 1 # Objective type
m.options.NODES = 2 # Collocation nodes
m.options.MAX_TIME = 10
##################################################################
# Create plot
plt.figure()
plt.ion()
plt.show()
filter_tc1 = []
def movefilter(predata, new, n):
if len(predata) < n:
predata.append(new)
else:
predata.pop(0)
predata.append(new)
return np.average(predata)
# Main Loop
start_time = time.time()
prev_time = start_time
last_Q1 = Q1_ss
try:
for i in range(1,loops):
# Sleep time
sleep_max = 2.0
print(time.time() - prev_time)
sleep = sleep_max - (time.time() - prev_time)
if sleep>=0.01:
time.sleep(sleep)
else:
time.sleep(0.01)
# Record time and change in time
t = time.time()
dt = t - prev_time
prev_time = t
tm[i] = t - start_time
# Read temperatures in Kelvin
curr_T1 = a.T1
last_T1 = curr_T1
avg_T1 = movefilter(filter_tc1, last_T1, 3)
T1[i] = curr_T1
###############################
### MPC CONTROLLER ###
###############################
cv.MEAS = avg_T1
# input setpoint with deadband +/- DT
DT = 0.5
cv.SPHI = Tsp1[i] + DT
cv.SPLO = Tsp1[i] - DT
# solve MPC
m.solve(disp=False)
# test for successful solution
if (m.options.APPSTATUS==1):
# retrieve the first Q value
Q1s[i] = mv.NEWVAL
print('Q1.NEWVAL', mv.NEWVAL)
with open(m.path+'//results.json') as f:
results = json.load(f)
else:
# not successful, set heater to zero
Q1s[i] = last_Q1
print('last_Q1', mv.NEWVAL)
print(m.path)
# Write output (0-100)
a.Q1(Q1s[i])
# Plot
plt.clf()
ax=plt.subplot(2,1,1)
ax.grid()
plt.plot(tm[0:i],T1[0:i],'ro',MarkerSize=3,label=r'$T_1$')
plt.plot(tm[0:i],Tsp1[0:i],'b-',MarkerSize=3,label=r'$T_1 Setpoint$')
cv_name = cv.NAME + '.bcv'
print(cv_name)
plt.plot(tm[i]+m.time,results[cv_name],'k-.',\
label=r'$T_1$ predicted',lw=3)
# plt.plot(tm[i]+m.time,results['tc1.tr_hi'],'k--',\
# label=r'$T_1$ trajectory')
# plt.plot(tm[i]+m.time,results['tc1.tr_lo'],'k--')
plt.ylabel('Temperature (degC)')
plt.legend(loc='best')
ax=plt.subplot(2,1,2)
ax.grid()
plt.plot(tm[0:i],Q1s[0:i],'r-',lw=3,label=r'$Q_1$')
plt.plot(tm[i]+m.time,mv.value,'k-.',\
label=r'$Q_1$ plan',lw=3)
# plt.plot(tm[0:i],Q2s[0:i],'b:',lw=3,label=r'$Q_2$')
plt.ylabel('Heaters')
plt.xlabel('Time (sec)')
plt.legend(loc='best')
plt.draw()
plt.pause(0.05)
# Turn off heaters
a.Q1(0)
a.Q2(0)
print('Shutting down')
a.close()
# Allow user to end loop with Ctrl-C
except KeyboardInterrupt:
# Disconnect from Arduino
a.Q1(0)
a.Q2(0)
print('Shutting down')
a.close()
# Make sure serial connection still closes when there's an error
except:
# Disconnect from Arduino
a.Q1(0)
a.Q2(0)
print('Error: Shutting down')
a.close()
- imaplib "command SELECT illegal in state NONAUTH" error
- Django - Login with Email
- Restart python thread with os.execl method
- llama-cpp-python not using NVIDIA GPU CUDA
- How to print multiple lines of text with Python
- List all contiguous sub-arrays
- How to list all existing loggers using python.logging module
- TRIO Lib queue get and put
- PyQt6 issue in fetching geometry of the window
- How to make a calculator function that uses the operator as an input?
- Using python range function for multiplication
- Notification after shut down on windows using python
- Python: Get unbound class method
- how to see if a word starts with char g-z(python 3)
- Error in anyjson setup command: use_2to3 is invalid
- What is gettext_lazy on django for?
- How do I create a multiline plot using seaborn?
- How to remove Empty Rows from a list in Python 3.5
- No module named 'PIL'
- Multiple inheritance metaclass conflict
- Python 3: Move email to trash by uid (imaplib)
- Python Module Not Found Error only for Compiled Script with Pyogrio Import
- How to add a state to a Selection field and show it ordered in the statusbar in Odoo 11?
- Breaking long method chains into multiple lines in Python
- Macro VS Micro VS Weighted VS Samples F1 Score
- Receiving Import Error: No Module named ***, but has __init__.py
- Hex string to signed int in Python
- How to convert a hex string to signed integer in Python
- server in a thread (Python3.9.0+aiohttp) : RuntimeError: can't register atexit after shutdown
- How to display a mathematical equation in a pretty way in Python 3