Fit Second Order with Optimization
 |  |  |  |  |
|---|
Fit parameters `K_p` and `\tau_p` from a first order process.
$$G_1(s)=\frac{K_p}{\tau_p \, s + 1}$$
The first order process is in parallel with another first order process.
$$G_2(s)=\frac{-2}{s+1}$$
The combined system has a second order response with data points sampled at 0.5 second intervals from a step input change of 2.0.
time, output (y) 0 , 0 0.5 , -0.2467 1 , -0.1677 1.5 , 0.0583 2 , 0.3341 2.5 , 0.6093 3 , 0.8604 3.5 , 1.0781 4 , 1.2613 4.5 , 1.412 5 , 1.5344 5.5 , 1.6328 6 , 1.7112 6.5 , 1.7734 7 , 1.8225 7.5 , 1.8611 8 , 1.8914 8.5 , 1.9152 9 , 1.9338 9.5 , 1.9484 10 , 1.9598
Use optimization to minimize the difference between a predicted response and the 21 measured values. Python and Excel have optimization solvers that may be useful for this exercise.
Scipy Optimize Minimize
import pandas as pd
import matplotlib.pyplot as plt
from scipy.optimize import minimize
# Import data file
# Column 1 = time (t)
# Column 2 = output (ymeas)
url = 'https://apmonitor.com/pdc/uploads/Main/data_2nd_order.txt'
data = pd.read_csv(url)
t = np.array(data['time'])
ymeas = np.array(data['output (y)'])
def model(p):
Kp = p[0]
taup = p[1]
# predicted values
ypred = 2.0 * Kp * (1.0-np.exp(-t/taup)) \
- 4.0 * (1-np.exp(-t))
return ypred
def objective(p):
ypred = model(p)
sse = sum((ymeas-ypred)**2)
return sse
# initial guesses for Kp and taup
p0 = [1.0,1.0]
# show initial objective
print('Initial SSE Objective: ' + str(objective(p0)))
# optimize Kp, taup
solution = minimize(objective,p0)
p = solution.x
# show final objective
print('Final SSE Objective: ' + str(objective(p)))
print('Kp: ' + str(p[0]))
print('taup: ' + str(p[1]))
# calculate new ypred
ypred = model(p)
# plot results
plt.figure()
plt.plot(t,ypred,'r-',linewidth=3,label='Predicted')
plt.plot(t,ymeas,'ko',linewidth=2,label='Measured')
plt.ylabel('Output')
plt.xlabel('Time')
plt.legend(loc='best')
plt.savefig('optimized.png')
plt.show()
GEKKO Optimization Suite
This source code takes a different approach than the Scipy.optimize.minimize or Excel solution where the explicit solution is matched to the data. In this case, the differential equations and the objective function are solved simultaneously. The two first order transfer functions are converted back to differential equation form.
$$\frac{Y_1(s)}{U(s)}=\frac{Kp}{\tau_p s + 1} \rightarrow \tau_p \frac{dy_1}{dt} + y_1 = K_p u$$
$$\frac{Y_2(s)}{U(s)}=\frac{-2}{\tau_p s + 1} \rightarrow \frac{dy_2}{dt} + y_2 = -2 u$$
The two responses `y_1` and `y_2` are computed separately and added together to predict the output `y`.
$$y=y_1+y_2$$
The sum of squared error objective function minimizes the difference between the predicted `y` and the measured `y_{meas}`.
$$\min_{K_p,\tau_p} \left(y-y_{meas}\right)^2$$
import pandas as pd
import matplotlib.pyplot as plt
from gekko import GEKKO
# Import data file
# Column 1 = time (t)
# Column 2 = output (ymeas)
url = 'https://apmonitor.com/pdc/uploads/Main/data_2nd_order.txt'
data = pd.read_csv(url)
t = np.array(data['time'])
ymeas = np.array(data['output (y)'])
# Create GEKKO model
m = GEKKO()
m.time = t
# Inputs
u = 2 # step input
ym = m.Param(ymeas)
# estimate parameters Kp and taup
Kp = m.FV()
taup = m.FV(lb=0)
# turn on parameters (STATUS=1) to let optimize change them
Kp.STATUS = 1
taup.STATUS = 1
# Equation 1
y1 = m.Var(0.0)
m.Equation(taup*y1.dt()==-y1+Kp*u)
# Equation 2
y2 = m.Var(0.0)
m.Equation(y2.dt()+y2==-2*u)
# Summation
y = m.Var(0.0)
m.Equation(y==y1+y2)
# Minimize Objective
m.Obj((y-ym)**2)
# Optimize Kp, taup
m.options.IMODE=5
m.options.NODES=3
m.solve()
# show final objective
print('Final SSE Objective: ' + str(m.options.OBJFCNVAL))
print('Kp: ' + str(Kp.value[0]))
print('taup: ' + str(taup.value[0]))
# plot results
plt.figure()
plt.plot(t,y.value,'r-',linewidth=3,label='Predicted')
plt.plot(t,ymeas,'ko',linewidth=2,label='Measured')
plt.ylabel('Output')
plt.xlabel('Time')
plt.legend(loc='best')
plt.savefig('optimized.png')
plt.show()
