Is there any example of data-driven model control in the course?
197 views
Skip to first unread message
Lifan Guo
Sep 11, 2018, 8:29:20 AM9/11/18
to apmonitor
Dear Prof.Hedengren,
I watched almost all of videos on your websites. It is awesome.
But, I wonder if there is any example of data-driven model with control?
It seems that all lectures include first principle model already, I am interested in how we can do system identification from data and do model predictive control on that.
I searched with google and I found some lectures from Professor Brunton, but there is not too much examples in python.
Many Thanks.
Regards.
John Hedengren
Sep 11, 2018, 9:46:37 AM9/11/18
to APM Google Groups
Dear Lifan,
There is an apm_id function for both APM MATLAB and APM Python that facilitates identification of ARX and Output Error models. We haven't added it yet for GEKKO (Python). Here is an example of system identification for MATLAB: https://github.com/APMonitor/apm_matlab/tree/master/example_id and an example for Python: https://github.com/APMonitor/apm_python/blob/master/example_lti_regression/lti_regression.py
For the course, there is a page for MIMO identification here: http://apmonitor.com/do/index.php/Main/ModelIdentification and the Temperature Control Lab has a few exercises for structured versus unstructured modeling: http://apmonitor.com/do/index.php/Main/AdvancedTemperatureControl
For empirical regression (not time series models), there are several other resources:
- Piece-wise Linear: http://apmonitor.com/wiki/index.php/Apps/PiecewiseLinear
- Nonlinear Regression: http://apmonitor.com/che263/index.php/Main/PythonDataRegression
- Another Nonlinear Regression: http://apmonitor.com/me575/index.php/Main/NonlinearRegression
- Regression examples in GEKKO including Neural Network:
https://apmonitor.com/wiki/index.php/Main/GekkoPythonOptimization
Let me know if you need additional resources.
Best regards,
John Hedengren
--
--
APMonitor user's group e-mail list.
- To post a message, send email to apmo...@googlegroups.com
- To unsubscribe, send email to apmonitor+unsubscribe@googlegroups.com
- Visit this group at http://groups.google.com/group/apmonitor
---
You received this message because you are subscribed to the Google Groups "apmonitor" group.
To unsubscribe from this group and stop receiving emails from it, send an email to apmonitor+unsubscribe@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
--
Best regards,
Lifan Guo
Sep 12, 2018, 6:31:55 AM9/12/18
to apmo...@googlegroups.com
Thank you~ It is really insightful.
I followed the python code in
with a small modification, since I don't have hardware.
I got a different result as below:
I wonder if this is the right result, cause it looks different with the one in the video.
My code is :
import tclab
TCLab = tclab.setup(connected=False, speedup=30)
with TCLab() as a:
print("Temperature 1: {0:0.2f} °C".format(a.T1))
print("Temperature 2: {0:0.2f} °C".format(a.T2))
# the rest part is as the same as the one online.
import tclab # pip install tclab
import numpy as np
import time
import matplotlib.pyplot as plt
from scipy.integrate import odeint
%matplotlib inline
#coding:utf-8
plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文标签
plt.rcParams['axes.unicode_minus']=False #用来正常显示负号
# FOPDT model
Kp = 0.5 # degC/%
tauP = 120.0 # seconds
thetaP = 10 # seconds (integer)
Tss = 23 # degC (ambient temperature)
Qss = 0 # % heater
# define energy balance model
def heat(x,t,Q):
# Parameters
Ta = 23 + 273.15 # K
U = 10.0 # W/m^2-K
m = 4.0/1000.0 # kg
Cp = 0.5 * 1000.0 # J/kg-K
A = 12.0 / 100.0**2 # Area in m^2
alpha = 0.01 # W / % heater
eps = 0.9 # Emissivity
sigma = 5.67e-8 # Stefan-Boltzman
# Temperature State
T = x[0]
# Nonlinear Energy Balance
dTdt = (1.0/(m*Cp))*(U*A*(Ta-T) \
+ eps * sigma * A * (Ta**4 - T**4) \
+ alpha*Q)
return dTdt
# save txt file
def save_txt(t,u1,u2,y1,y2,sp1,sp2):
data = np.vstack((t,u1,u2,y1,y2,sp1,sp2)) # vertical stack
data = data.T # transpose data
top = 'Time (sec), Heater 1 (%), Heater 2 (%), ' \
+ 'Temperature 1 (degC), Temperature 2 (degC), ' \
+ 'Set Point 1 (degC), Set Point 2 (degC)'
np.savetxt('data.txt',data,delimiter=',',header=top,comments='')
print('LED On')
a.LED(100)
# Run time in minutes
run_time = 10.0
# Number of cycles
loops = int(60.0*run_time)
tm = np.zeros(loops)
# Temperature (K)
Tsp1 = np.ones(loops) * 23.0 # set point (degC)
T1 = np.ones(loops) * a.T1 # measured T (degC)
Tsp2 = np.ones(loops) * 23.0 # set point (degC)
T2 = np.ones(loops) * a.T2 # measured T (degC)
# Predictions
Tp = np.ones(loops) * a.T1
error_eb = np.zeros(loops)
Tpl = np.ones(loops) * a.T1
error_fopdt = np.zeros(loops)
# impulse tests (0 - 100%)
Q1 = np.ones(loops) * 0.0
Q2 = np.ones(loops) * 0.0
Q1[10:110] = 50.0 # step up for 100 sec
Q1[200:300] = 90.0 # step up for 100 sec
Q1[400:500] = 70.0 # step up for 100 sec
print('Running Main Loop. Ctrl-C to end.')
print(' Time Q1 Q2 T1 T2')
print('{:6.1f} {:6.2f} {:6.2f} {:6.2f} {:6.2f}'.format(tm[0], \
Q1[0], \
Q2[0], \
T1[0], \
T2[0]))
# Create plot
plt.figure(figsize=(10,7))
plt.ion()
plt.show()
# Main Loop
start_time = time.time()
prev_time = start_time
try:
for i in range(1,loops):
# Sleep time
sleep_max = 1.0
sleep = sleep_max - (time.time() - prev_time)
if sleep>=0.01:
time.sleep(sleep-0.01)
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
T1[i] = a.T1
T2[i] = a.T2
# Simulate one time step with Energy Balance
Tnext = odeint(heat,Tp[i-1]+273.15,[0,dt],args=(Q1[i-1],))
Tp[i] = Tnext[1]-273.15
error_eb[i] = error_eb[i-1] + abs(Tp[i]-T1[i])
# Simulate one time step with FOPDT model
z = np.exp(-dt/tauP)
Tpl[i] = (Tpl[i-1]-Tss) * z \
+ (Q1[max(0,i-int(thetaP)-1)]-Qss)*(1-z)*Kp \
+ Tss
error_fopdt[i] = error_fopdt[i-1] + abs(Tpl[i]-T1[i])
# Write output (0-100)
a.Q1(Q1[i])
a.Q2(Q2[i])
# Print line of data
print('{:6.1f} {:6.2f} {:6.2f} {:6.2f} {:6.2f}'.format(tm[i], \
Q1[i], \
Q2[i], \
T1[i], \
T2[i]))
# Plot
plt.clf()
ax=plt.subplot(3,1,1)
ax.grid()
plt.plot(tm[0:i],T1[0:i],'ro',label=r'$T_1$ measured')
plt.plot(tm[0:i],Tp[0:i],'k-',label=r'$T_1$ energy balance')
plt.plot(tm[0:i],Tpl[0:i],'g:',label=r'$T_1$ FOPDT')
plt.plot(tm[0:i],T2[0:i],'bx',label=r'$T_2$ measured')
plt.ylabel('Temperature (degC)')
plt.legend(loc=2)
ax=plt.subplot(3,1,2)
ax.grid()
plt.plot(tm[0:i],error_eb[0:i],'k-',label='Energy Balance')
plt.plot(tm[0:i],error_fopdt[0:i],'g:',label='Linear')
plt.ylabel('Cumulative Error')
plt.legend(loc='best')
ax=plt.subplot(3,1,3)
ax.grid()
plt.plot(tm[0:i],Q1[0:i],'r-',label=r'$Q_1$')
plt.plot(tm[0:i],Q2[0:i],'b:',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)
# Save text file
save_txt(tm[0:i],Q1[0:i],Q2[0:i],T1[0:i],T2[0:i],Tsp1[0:i],Tsp2[0:i])
# Save figure
plt.savefig('test_Models.png')
except KeyboardInterrupt:
# Disconnect from Arduino
a.Q1(0)
a.Q2(0)
print('Shutting down')
a.close()
save_txt(tm[0:i],Q1[0:i],Q2[0:i],T1[0:i],T2[0:i],Tsp1[0:i],Tsp2[0:i])
plt.savefig('test_Models.png')
# 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()
save_txt(tm[0:i],Q1[0:i],Q2[0:i],T1[0:i],T2[0:i],Tsp1[0:i],Tsp2[0:i])
plt.savefig('test_Models.png')
raise
Many thanks.
On Tue, Sep 11, 2018 at 9:46 PM John Hedengren <jo...@apmonitor.com> wrote:
Dear Lifan,There is an apm_id function for both APM MATLAB and APM Python that facilitates identification of ARX and Output Error models. We haven't added it yet for GEKKO (Python). Here is an example of system identification for MATLAB: https://github.com/APMonitor/apm_matlab/tree/master/example_id and an example for Python: https://github.com/APMonitor/apm_python/blob/master/example_lti_regression/lti_regression.pyFor the course, there is a page for MIMO identification here: http://apmonitor.com/do/index.php/Main/ModelIdentification and the Temperature Control Lab has a few exercises for structured versus unstructured modeling: http://apmonitor.com/do/index.php/Main/AdvancedTemperatureControlFor empirical regression (not time series models), there are several other resources:
- Piece-wise Linear: http://apmonitor.com/wiki/index.php/Apps/PiecewiseLinear
- Nonlinear Regression: http://apmonitor.com/che263/index.php/Main/PythonDataRegression
- Another Nonlinear Regression: http://apmonitor.com/me575/index.php/Main/NonlinearRegression
- Regression examples in GEKKO including Neural Network: https://apmonitor.com/wiki/index.php/Main/GekkoPythonOptimization
Let me know if you need additional resources.Best regards,John Hedengren
On Mon, Sep 10, 2018 at 10:06 PM, Lifan Guo <lifa...@gmail.com> wrote:
Dear Prof.Hedengren,I watched almost all of videos on your websites. It is awesome.But, I wonder if there is any example of data-driven model with control?It seems that all lectures include first principle model already, I am interested in how we can do system identification from data and do model predictive control on that.I searched with google and I found some lectures from Professor Brunton, but there is not too much examples in python.Many Thanks.Regards.
--
--
APMonitor user's group e-mail list.
- To post a message, send email to apmo...@googlegroups.com
- To unsubscribe, send email to apmonitor+...@googlegroups.com
- Visit this group at http://groups.google.com/group/apmonitor
---
You received this message because you are subscribed to the Google Groups "apmonitor" group.
To unsubscribe from this group and stop receiving emails from it, send an email to apmonitor+...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
--
--
APMonitor user's group e-mail list.
- To post a message, send email to apmo...@googlegroups.com
- To unsubscribe, send email to apmonitor+...@googlegroups.com
- Visit this group at http://groups.google.com/group/apmonitor
---
You received this message because you are subscribed to the Google Groups "apmonitor" group.
To unsubscribe from this group and stop receiving emails from it, send an email to apmonitor+...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
Lifan Guo
Sep 12, 2018, 6:31:55 AM9/12/18
to apmo...@googlegroups.com
oh, I think the speedup part make the result difference.
John Hedengren
Sep 12, 2018, 6:42:35 PM9/12/18
to APM Google Groups
Lifan,
Yes, I think you'll need to put in an artificial pause in order for the simulated and measured values to match. The TCLabs are now available on Amazon (https://www.amazon.com/gp/product/B07GMFWMRY) as well as through the Paypal link (https://apmonitor.com/pdc/index.php/Main/PurchaseLabKit) if you'd like to get a physical device. You should see a time constant of about 60-100 sec.
Best regards,
John Hedengren
- To unsubscribe, send email to apmonitor+unsubscribe@googlegroups.com
- Visit this group at http://groups.google.com/group/apmonitor
---
You received this message because you are subscribed to the Google Groups "apmonitor" group.
To unsubscribe from this group and stop receiving emails from it, send an email to apmonitor+unsubscribe@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
--
--
APMonitor user's group e-mail list.
- To post a message, send email to apmo...@googlegroups.com
- To unsubscribe, send email to apmonitor+unsubscribe@googlegroups.com
- Visit this group at http://groups.google.com/group/apmonitor
---
You received this message because you are subscribed to the Google Groups "apmonitor" group.
To unsubscribe from this group and stop receiving emails from it, send an email to apmonitor+unsubscribe@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
--
--
- To unsubscribe, send email to apmonitor+unsubscribe@googlegroups.com
- Visit this group at http://groups.google.com/group/apmonitor
---
You received this message because you are subscribed to the Google Groups "apmonitor" group.
To unsubscribe from this group and stop receiving emails from it, send an email to apmonitor+unsubscribe@googlegroups.com.
Lifan Guo
Sep 17, 2018, 8:08:38 AM9/17/18
to apmo...@googlegroups.com
Dear Prof.Hedengren,
I am working on understanding Moving Horizon Estimation part from the github code : https://github.com/APMonitor/arduino/tree/master/5_Moving_Horizon_Estimation/1st_order_linear/Python
I can run it successfully, but I have some questions regarding with the code:
1, where could I find the detail of how mhe function finds Kp,tau, TC_ss,TC given T_meas and Q1? It looks like a black box for me to understand
2, I find the equations in model.apm file and I think it estimated by first oder system to estimate the parameters. So I am trying use the parameters
I learned from MHE function to get the same value. But, it fails.
I am attaching the code and result here:
the function I use is the one in "https://github.com/APMonitor/arduino/tree/master/5_Moving_Horizon_Estimation/1st_order_linear/Python"
then I get estimated KP, TAU, TC_ss and TC.
Then I am trying to use this value and first oder system in python to get the same value
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
from scipy.integrate import odeint
# Simulate taup * dy/dt = -y + K*u
Kp = 0.45
taup = 158
# (3) ODE Integrator
def model3(y,t):
u = 100
return (-y + Kp * u)/taup
t3 = np.arange(0,120,2)
y3 = odeint(model3,0,t3)
print(y3)
I wonder if this is correct? cause I get y3 [-1]=23.67613127 while MHE result = 24.092
Best
- To unsubscribe, send email to apmonitor+...@googlegroups.com
- Visit this group at http://groups.google.com/group/apmonitor
---
You received this message because you are subscribed to the Google Groups "apmonitor" group.
To unsubscribe from this group and stop receiving emails from it, send an email to apmonitor+...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
--
--
APMonitor user's group e-mail list.
- To post a message, send email to apmo...@googlegroups.com
- To unsubscribe, send email to apmonitor+...@googlegroups.com
- Visit this group at http://groups.google.com/group/apmonitor
---
You received this message because you are subscribed to the Google Groups "apmonitor" group.
To unsubscribe from this group and stop receiving emails from it, send an email to apmonitor+...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
--
----
--
APMonitor user's group e-mail list.
- To post a message, send email to apmo...@googlegroups.com
- To unsubscribe, send email to apmonitor+...@googlegroups.com
- Visit this group at http://groups.google.com/group/apmonitor
---
You received this message because you are subscribed to the Google Groups "apmonitor" group.
To unsubscribe from this group and stop receiving emails from it, send an email to apmonitor+...@googlegroups.com.
--
--
APMonitor user's group e-mail list.
- To post a message, send email to apmo...@googlegroups.com
- To unsubscribe, send email to apmonitor+...@googlegroups.com
- Visit this group at http://groups.google.com/group/apmonitor
---
You received this message because you are subscribed to the Google Groups "apmonitor" group.
To unsubscribe from this group and stop receiving emails from it, send an email to apmonitor+...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
John Hedengren
Sep 17, 2018, 11:18:08 AM9/17/18
to apmo...@googlegroups.com
Lifan,
MHE updates the parameters each cycle as new measurements arrive. There is more detail on the MHE algorithm at the following links:
MHE Estimator Objective
http://apmonitor.com/do/index.php/Main/EstimatorObjective
Orthogonal Collocation (Direct Transcription) of DAE -> NLP
http://apmonitor.com/do/index.php/Main/OrthogonalCollocation
More Details on MHE
http://apmonitor.com/wiki/uploads/Main/APMonitor_2014.pdf
If you are interested in additional details, there is the Dynamic Optimization course that will be starting again in Jan 2019.
Best regards,
John Hedengren
Lifan Guo
Sep 19, 2018, 7:29:40 AM9/19/18
to apmo...@googlegroups.com
Dear Professor,
I read all links you provided above, and I have one question in the example shown in the link: "http://apmonitor.com/do/index.php/Main/EstimatorObjective"
I read the code in Motion Estimation from GPS Position (1D Example). I am trying to use scipy to solve the ode function in this example and
try to use scipy to do optimization based on square error. ( as I want to get better understanding by comparing scipy solution and apmonitor solution). However,
I have difficulty to solve ode function in this example:
dX1(t) /dt = X2(t) dX2(t) /dt = u(t)
I assign a random number to U. I know it is absolutely wrong, but could you please show me the correct code to do ODE based on the function?
my code is:
def dP_dt(P,t,u):
return [P[1], u]
ts = np.linspace(0, 12, 100)
u=np.linspace(0, 12, 100)
P0 = [0, 1.0]
Ps = odeint(dP_dt, P0, ts,args=(1)) # absolutely wrong,but how to correct this?
Best Regards.
John Hedengren
Sep 19, 2018, 3:10:42 PM9/19/18
to apmo...@googlegroups.com
Lifan,
Here are some resources:
· Solve ODEs with ODEINT: https://apmonitor.com/pdc/index.php/Main/SolveDifferentialEquations
· Solve ODEs with GEKKO: https://apmonitor.com/pdc/index.php/Main/PythonDifferentialEquations
I don’t have a sequential method example of MHE like you are trying to develop with an ODE integrator but there is an MPC sequential method here (developed by Junho Park):
http://apmonitor.com/pdc/index.php/Main/ModelPredictiveControl
You may be able to start with this as a template. In general, a sequential method is much slower than a simultaneous method for most systems. It is a nice way to learn the methods, however. If you want to implement a simultaneous method, there is information here: http://apmonitor.com/do/index.php/Main/OrthogonalCollocation
Reply all
Reply to author
Forward
0 new messages