System of FOPDT equations in GEKKO - use more than one input

Viewed 309

I want to use two inputs or more to create a more precise estimation of a variable. I already estimated it using only one input and one FOPDT equation, but when I try to add one more input and the respective k, tau and theta, along with another equation, i get "Solution Not Found" error. Can I create a system of equations this way?

More details about the solver output below. Even though I added more variables than equations to use more than one input, this made my problem have negative degrees of freedom.

--------- APM Model Size ------------
 Each time step contains
   Objects      :            2
   Constants    :            0
   Variables    :           15
   Intermediates:            0
   Connections  :            4
   Equations    :            6
   Residuals    :            6
 
 Number of state variables:           1918
 Number of total equations: -         2151
 Number of slack variables: -            0
 ---------------------------------------
 Degrees of freedom       :           -233
 
 * Warning: DOF <= 0
 **********************************************
 Dynamic Estimation with Interior Point Solver
 **********************************************
  
  
 Info: Exact Hessian

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt
******************************************************************************

This is Ipopt version 3.12.10, running with linear solver ma57.

Number of nonzeros in equality constraint Jacobian...:     6451
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:     1673

Exception of type: TOO_FEW_DOF in file "IpIpoptApplication.cpp" at line 891:
 Exception message: status != TOO_FEW_DEGREES_OF_FREEDOM evaluated false: Too few degrees of freedom (rethrown)!

EXIT: Problem has too few degrees of freedom.
 
 An error occured.
 The error code is          -10
 
 
 ---------------------------------------------------
 Solver         :  IPOPT (v3.12)
 Solution time  :   2.279999999154825E-002 sec
 Objective      :   0.000000000000000E+000
 Unsuccessful with error code            0
 ---------------------------------------------------
 
 Creating file: infeasibilities.txt
 Use command apm_get(server,app,'infeasibilities.txt') to retrieve file
 @error: Solution Not Found

And here's the code

from gekko import GEKKO
import numpy as np
import pandas as pd
import plotly.express as px 

d19jc = d19jc.dropna()

d19jcSlice = d19jc.loc['2019-10-22 05:30:00':'2019-10-22 09:30:00'] #jc22

d19jcSlice.index = pd.to_datetime(d19jcSlice.index)
d19jcSliceGroupMin = d19jcSlice.groupby(pd.Grouper(freq='T')).mean()
data = d19jcSliceGroupMin.dropna()

xdf1 = data['Cond_PID_SP']
xdf2 = data['Front_PID_SP']
ydf1 = data['Cond_Center_Top_TC']

xms1 = pd.Series(xdf1)
xm1 = np.array(xms1)
xms2 = pd.Series(xdf2)
xm2 = np.array(xms2)
yms = pd.Series(ydf1)
ym = np.array(yms)

xm_r = len(xm1)
tm = np.linspace(0,xm_r-1,xm_r)

m = GEKKO()
m.options.IMODE=5
m.time = tm; time = m.Var(0); m.Equation(time.dt()==1)

k1 = m.FV(lb=0.1,ub=5); k1.STATUS=1
tau1 = m.FV(lb=1,ub=300); tau1.STATUS=1
theta1 = m.FV(lb=0,ub=30); theta1.STATUS=1
k2 = m.FV(lb=0.1,ub=5); k2.STATUS=1
tau2 = m.FV(lb=1,ub=300); tau2.STATUS=1
theta2 = m.FV(lb=0,ub=30); theta2.STATUS=1

# create cubic spline with t versus u
uc1 = m.Var(xm1); tc1 = m.Var(tm); m.Equation(tc1==time-theta1)
m.cspline(tc1,uc1,tm,xm1,bound_x=False)
# create cubic spline with t versus u
uc2 = m.Var(xm2); tc2 = m.Var(tm); m.Equation(tc2==time-theta2)
m.cspline(tc2,uc2,tm,xm2,bound_x=False)

x1 = m.Param(value=xm1)
x2 = m.Param(value=xm2)
y = m.Var(value=ym)
yObj = m.Param(value=ym)

m.Equation(tau1*y.dt()+(y-ym[0])==k1 * (uc1-xm1[0]))
m.Equation(tau2*y.dt()+(y-ym[0])==k2 * (uc2-xm2[0]))
m.Minimize((y-yObj)**2)
m.options.EV_TYPE=2

print('solve start')
m.solve(disp=True)

print('k1: ', k1.value[0])
print('tau1: ',  tau1.value[0])
print('theta1: ', theta1.value[0])
print('k2: ', k2.value[0])
print('tau2: ',  tau2.value[0])
print('theta2: ', theta2.value[0])

df_plot = pd.DataFrame({'DateTime' : data.index,
                        'Cond_Center_Top_TC' : np.array(yObj),
                        'Fit Cond_Center_Top_TC - Train' : np.array(y),
figGekko = px.line(df_plot, 
                   x='DateTime',
                   y=['Cond_Center_Top_TC','Fit Cond_Center_Top_TC - Train'],
                   labels={"value": "Degrees Celsius"},
                   title = "(Cond_PID_SP & Front_PID_SP) -> Cond_Center_Top_TC JC only - Train")
figGekko.update_layout(legend_title_text='')
figGekko.show()
2 Answers

You currently have only one variable and two equations.

y = m.Var(value=ym)
yObj = m.Param(value=ym)

m.Equation(tau1*y.dt()+(y-ym[0])==k1 * (uc1-xm1[0]))
m.Equation(tau2*y.dt()+(y-ym[0])==k2 * (uc2-xm2[0]))

You need to create two separate variables y1 and y2 for the two equations.

y1 = m.Var(value=ym)
y2 = m.Var(value=ym)
yObj = m.Param(value=ym)

m.Equation(tau1*y1.dt()+(y1-ym[0])==k1 * (uc1-xm1[0]))
m.Equation(tau2*y2.dt()+(y2-ym[0])==k2 * (uc2-xm2[0]))

You may also need to create two separate measurement vectors for ym1 and ym2 if you have different data for each.

Edit: Multiple Input, Single Output (MISO) System

For a MISO system, you need to adjust the equation as shown in the Process Dynamics and Control course.

y = m.Var(value=ym)
yObj = m.Param(value=ym)

m.Equation(tau*y.dt()+(y-ym[0])==k1 * (uc1-xm1[0]) + k2 * (uc2-xm2[0]))

There is only one time constant in this form. If they need different time constants, you can also add the two outputs together with:

y1 = m.Var(value=ym[0])
y2 = m.Var(value=ym[0])
y  = m.Var(value=ym)
yObj = m.Param(value=ym)

m.Equation(tau1*y1.dt()+(y1-ym[0])==k1 * (uc1-xm1[0]))
m.Equation(tau2*y2.dt()+(y2-ym[0])==k2 * (uc2-xm2[0]))
m.Equation(y==y1+y2)

In this case, y is the variable that has data and y1 and y2 are unmeasured states.

Please see the below example code for a FOPDT model that has two inputs and one output. Both the transfer function models have different FOPDT parameters, including different deadtimes as well.
It is still in a dynamic simulation mode (IMODE=4), but you can start from this modifying a little bit toward the dynamic estimation mode (IMODE=5) and MPC mode (IMODE=6) later.

from gekko import GEKKO
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

tf = 100 
npt = 101 
t = np.linspace(0,tf,npt)
u1 = np.zeros(npt)
u2 = np.zeros(npt)
u1[10:] = 5
u2[40:] = -5

m = GEKKO(remote=True)
m.time = t 
time = m.Var(0) 
m.Equation(time.dt()==1)

K1 = m.FV(1,lb=0,ub=1);      K1.STATUS=1
tau1 = m.FV(5, lb=1,ub=300);  tau1.STATUS=1
theta1 = m.FV(10, lb=2,ub=30); theta1.STATUS=1

K2 = m.FV(0.5,lb=0,ub=1);      K2.STATUS=1
tau2 = m.FV(10, lb=1,ub=300);  tau2.STATUS=1
theta2 = m.FV(20, lb=2,ub=30); theta2.STATUS=1

uc1 = m.Var(u1) 
uc2 = m.Var(u2) 
tc1 = m.Var(t) 
tc2 = m.Var(t)
m.Equation(tc1==time-theta1)
m.Equation(tc2==time-theta2)
m.cspline(tc1,uc1,t,u1,bound_x=False)
m.cspline(tc2,uc2,t,u2,bound_x=False)

yp = m.Var() 
yp1 = m.Var()
yp2 = m.Var()
m.Equation(yp1.dt() == -yp1/tau1 + K1*uc1/tau1) 
m.Equation(yp2.dt() == -yp2/tau2 + K2*uc2/tau2)
m.Equation(yp == yp1 + yp2)

m.options.IMODE=4
m.solve()

print('K1: ', K1.value[0])
print('tau1: ',  tau1.value[0])
print('theta1: ', theta1.value[0])
print('')

print('K2: ', K2.value[0])
print('tau2: ',  tau2.value[0])
print('theta2: ', theta2.value[0])

plt.figure()
plt.subplot(2,1,1)
plt.plot(t,u1)
plt.plot(t,u2)
plt.legend([r'u1', r'u2'])
plt.ylabel('Inputs 1 & 2')
plt.subplot(2,1,2)
plt.plot(t,yp)
plt.legend([r'y1'])
plt.ylabel('Output')
plt.xlabel('Time')
plt.savefig('sysid.png')
plt.show()

enter image description here

K1:  1.0
tau1:  5.0
theta1:  10.0

K2:  0.5
tau2:  10.0
theta2:  20.0
Related