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 numpy as np
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 numpy as np
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()