pythongekkompc

System identification using an ARX model with Gekko


The following is related to this question: Why doesn't the Gekko solver adapt to variations in the system?

What I still don't understand is why, sometimes, when the outside temperature rises, the inside temperature remains constant. Normally, it should also increase since beta remains constant. Here are examples of this:

example 1

example 2

In the two images, we can see that on the second day, the outside temperature increases and β remains constant, but the inside temperature decreases.

I checked my data and found that it had been taken every 2 seconds, which represents 43,200 measurements per day. However, during the system identification with the ARX model, I had set the sampling time to 300 sec.

I ran the corrected MPC code with 2 sec sampling from my previous question.

# Import library
import numpy as np
import pandas as pd
import time
from gekko import GEKKO
from numpy import array
import matplotlib.pyplot as plt

# import data
tx = pd.read_csv(r'C:\Users\mpc\tx_b.csv')
tz = pd.read_csv(r'C:\Users\mpc\tz.csv')

data = pd.concat([tx,tz],axis=1)
data_1 = data[0:8596800]
data_2 = data[8596800:]

#%%  Initialize Model
m = GEKKO(remote=False)

# system identification
ts = 2
t = np.arange(0,len(data_1)*ts, ts)
u_id = data_1[['Tx_i','beta_i']]
y_id = data_1[['Tz_i']]

#meas : the time-series next step is predicted from prior measurements as in ARX

na=10; nb=10 # ARX coefficients
print('Identify model')
start = time.time()
yp,p,K = m.sysid(t,u_id,y_id,na,nb,objf=100,scale=False,diaglevel=0,pred='meas')
print('temps de prediction :'+str(time.time()-start)+'s')

#%% parametres 
K = array([[  0.93688819, -12.22410568]])
p = {'a': array([[ 1.08945931],
        [-0.00243571],
        [-0.00247112],
        [-0.00273341],
        [-0.00296342],
        [-0.00319516],
        [-0.00343794],
        [-0.00366398],
        [-0.00394255],
        [-0.06661506]]),
 'b': array([[[-0.05134201, -0.01035174],
         [ 0.00170311, -0.01551259],
         [ 0.00172715, -0.01178932],
         [ 0.00178147, -0.01051817],
         [ 0.00184694, -0.00821511],
         [ 0.00192371, -0.00570574],
         [ 0.00201409, -0.00344425],
         [ 0.00210016, -0.0014708 ],
         [ 0.00222189,  0.00021622],
         [ 0.03789636,  0.04235503]]]),
 'c': array([0.0266222])}


#%% I used the last day's external temperature data as a disturbance.
T_externel = data_2[["Tx_i"]].values

m = GEKKO(remote=False)
m.y = m.Array(m.CV,1)
m.u = m.Array(m.MV,2)
m.arx(p,m.y,m.u)

# rename CVs
m.T = m.y[0]

# rename MVs
m.beta = m.u[1]

# distrubance
m.d = m.u[0]

# distrubance and parametres
m.d = m.Param(T_externel[0])

m.bias = m.Param(0)
m.Tb = m.CV()
m.Equation(m.Tb==m.T+m.bias)

# steady state initialization
m.options.IMODE = 1
m.solve(disp=False)

# set up MPC
m.d.value = T_externel

m.options.IMODE   = 6 # MPC
m.options.CV_TYPE = 1 # the objective is an l1-norm (region)
m.options.NODES   = 2 # Collocation nodes
m.options.SOLVER  = 3 # IPOPT
m.time = np.arange(0,len(T_externel)*2,2) # step time = 300s

# Manipulated variables
m.beta.STATUS = 1  # calculated by the optimizer
m.beta.FSTATUS = 0 # use measured value
m.beta.DCOST = 0.0 # Delta cost penalty for MV movement
m.beta.UPPER = 1.0 # Upper bound
m.beta.LOWER = 0.0 # Lower bound
m.beta.MV_STEP_HOR = 1
m.beta.value = 0

# Controlled variables
m.Tb.STATUS = 1        # drive to set point
m.Tb.FSTATUS = 0       # receive measurement
m.Tb.SPHI = 17.5         # set point high level
m.Tb.SPLO = 16.5         # set point low level
m.Tb.WSPHI = 100       # set point high priority
m.Tb.WSPLO = 100       # set point low priority

T_MEAS = 20
m.Tb.value = T_MEAS
m.bias.value = T_MEAS - m.T.value[0]

m.options.SOLVER = 3
m.solve(disp=False)

if m.options.APPSTATUS == 1:
   # Retrieve new values
   beta  = m.beta.NEWVAL
   for i in range(43200):
       print(i)
else:
   # Solution failed
   beta  = 0.0

# Plot the results
plt.figure(figsize=(8,3.5))
plt.subplot(3,1,1)
plt.plot(m.time,m.Tb.value,'r-',label=r'$T_{int}$')
plt.plot([0,m.time[-1]],[m.Tb.SPHI,m.Tb.SPHI],'k--',label='Upper Bound')
plt.plot([0,m.time[-1]],[m.Tb.SPLO,m.Tb.SPLO],'k--',label='Lower Bound')
plt.legend(loc=1); plt.grid()
plt.ylabel('Tin (°C)')
plt.subplot(3,1,2)
plt.plot(m.time,m.d.value,'g:',label=r'$T_{ext}$')
plt.ylabel('Tex (°C)')
plt.subplot(3,1,3)
plt.step(m.time,m.beta.value,'b--',label=r'$\beta$')
plt.ylabel('optimal control')
plt.xlabel('Time (sec)')
plt.legend(loc=1); plt.grid()
plt.savefig('results7.png',dpi=300)
plt.show()

For the external temperature over one day, I took one day's worth of data (data_2) and concatenated it with itself to create a profile for two consecutive days.

Here are the plots for the two days

I haven't solved the problem yet.


Solution

  • There is a small error where m.d is redefined and the connection is broken with the m.arx() model. Here is a correction:

    # distrubance and parametres
    #m.d = m.Param(T_externel[0])
    m.d.value = T_externel[0]
    

    There is also a dynamics issue. The response to beta (and likely external temperature) is very slow as shown in the step response here:

    slow response

    This could be from using a different sampling time interval for the SysID (System Identification) and the controller. This is especially the case if you used a short time interval for system identification and use a longer time interval for the MPC application. Because it is a discrete ARX model, the dynamics will take longer for the inputs to affect the outputs.

    Response with Additional Information

    Thanks for confirming that the identification was performed with 2 second data and the MPC application with 300 second time intervals. This is the source of one of the issues. One method is to use 2 sec for both (as you have shown) or else use 300 sec for both. I recommend going to 300 sec intervals for both because of the slow dynamics and need to calculate over a long time horizon. There is a way to reduce the system identification data to 300 second intervals with a rolling average and sampling every 150 data points with davg.iloc[::150]. Here are the first 500 points of the identification.

    results

    # Import library
    import numpy as np
    import pandas as pd
    import time
    from gekko import GEKKO
    from numpy import array
    import matplotlib.pyplot as plt
    
    # import data
    tx = pd.read_csv('tx_b.csv')
    tz = pd.read_csv('tz.csv')
    data = pd.concat([tx,tz],axis=1)
    
    # reduce from 2 sec to 5 min intervals
    # calculate rolling average
    davg = data.rolling(window=150).mean()
    # sample every 150 cycles (2 sec to 5 min)
    dnew = davg.iloc[::150]
    dnew.dropna(inplace=True)
    
    slice = int(len(dnew)*0.99)
    data_1 = dnew.iloc[:slice].copy() # 99% for training
    data_2 = dnew.iloc[slice:].copy() # 1% for testing
    
    #%%  Initialize Model
    m = GEKKO(remote=False)
    
    # system identification
    ts = 300
    t = np.arange(0,len(data_1)*ts, ts)
    u_id = data_1[['Tx_i','beta_i']]
    y_id = data_1[['Tz_i']]
    
    #meas : the time-series next step is predicted from prior measurements as in ARX
    
    na=10; nb=10 # ARX coefficients
    print('Identify model')
    start = time.time()
    yp,p,K = m.sysid(t,u_id,y_id,na,nb,objf=100,scale=False,diaglevel=0,pred='meas')
    print('temps de prediction :'+str(time.time()-start)+'s')
    
    # new parameters
    print(f'Parameters: {p}')
    
    # new gains
    print(f'Gains: {K}')
    
    # Plot the results
    nview = 500
    plt.figure(figsize=(8,3.5))
    plt.subplot(3,1,1)
    plt.plot(t[0:nview],data_1[['Tz_i']].iloc[0:nview],'r-',label=r'$T_{meas}$')
    plt.plot(t[0:nview],yp[0:nview],'b--',label=r'$T_{pred}$')
    plt.legend(loc=1); plt.grid()
    plt.ylabel('Tin (°C)')
    plt.subplot(3,1,2)
    plt.plot(t[0:nview],data_1['Tx_i'].iloc[0:nview],'g:',label=r'$T_{ext}$')
    plt.ylabel('Tex (°C)')
    plt.legend(loc=1); plt.grid()
    plt.subplot(3,1,3)
    plt.step(t[0:nview],data_1['beta_i'].iloc[0:nview],'b--',label=r'$\beta$')
    plt.ylabel('optimal control')
    plt.xlabel('Time (sec)')
    plt.legend(loc=1); plt.grid()
    plt.savefig('results_id.png',dpi=300)
    plt.show()
    

    The new gains look good:

    Gains: [[  0.97529892 -14.28674666]]
    

    The new parameters are:

    Parameters:
    {'a': array([[ 0.65059382],
           [-0.03027903],
           [ 0.06080624],
           [-0.02233488],
           [ 0.06877417],
           [ 0.37139369],
           [-0.22535397],
           [-0.01737772],
           [-0.03430238],
           [ 0.02420281]]),
     'b': array([[[ 3.71568256e-01, -8.09407185e+00],
            [-7.62375332e-02,  4.97505860e+00],
            [-1.34657379e-01, -1.56110911e+00],
            [ 7.32255512e-02,  1.89452386e+00],
            [-2.50583697e-02, -1.68565102e+00],
            [-2.06554465e-01,  2.03970469e+00],
            [ 1.50625153e-01,  8.79190526e-01],
            [-9.20974935e-03, -1.34505749e-01],
            [ 1.58409954e-03, -6.32745151e-01],
            [ 4.79076587e-03,  1.21199721e-01]]]),
     'c': array([2.08698176])}
    

    Here is a step response with T_externel.

    step response

    # Import library
    import numpy as np
    import pandas as pd
    import time
    from gekko import GEKKO
    from numpy import array
    import matplotlib.pyplot as plt
    
    #%% parameters
    p = {'a': array([[ 0.65059382],
           [-0.03027903],
           [ 0.06080624],
           [-0.02233488],
           [ 0.06877417],
           [ 0.37139369],
           [-0.22535397],
           [-0.01737772],
           [-0.03430238],
           [ 0.02420281]]),
         'b': array([[[ 3.71568256e-01, -8.09407185e+00],
            [-7.62375332e-02,  4.97505860e+00],
            [-1.34657379e-01, -1.56110911e+00],
            [ 7.32255512e-02,  1.89452386e+00],
            [-2.50583697e-02, -1.68565102e+00],
            [-2.06554465e-01,  2.03970469e+00],
            [ 1.50625153e-01,  8.79190526e-01],
            [-9.20974935e-03, -1.34505749e-01],
            [ 1.58409954e-03, -6.32745151e-01],
            [ 4.79076587e-03,  1.21199721e-01]]]),
         'c': array([2.08698176])}
    
    m = GEKKO(remote=False)
    m.y = m.Array(m.CV,1)
    m.u = m.Array(m.MV,2)
    m.arx(p,m.y,m.u)
    
    # rename CVs
    m.T = m.y[0]
    
    # rename MVs
    m.beta = m.u[1]
    
    # distrubance
    m.d = m.u[0]
    
    # distrubance and parametres
    m.d.value = 0
    
    m.bias = m.Param(0)
    m.Tb = m.CV()
    m.Equation(m.Tb==m.T+m.bias)
    
    # steady state initialization
    m.options.IMODE = 1
    m.solve(disp=False)
    
    # set up MPC
    n = 120 # 10 hours at 5 min sampling rate
    T_externel = np.zeros(n)
    T_externel[10:] = 15
    
    m.d.value = T_externel
    
    m.options.IMODE   = 6 # MPC
    m.options.CV_TYPE = 1 # the objective is an l1-norm (region)
    m.options.NODES   = 2 # Collocation nodes
    m.options.SOLVER  = 3 # IPOPT
    m.time = np.arange(0,len(T_externel)*300,300) # step time = 300s
    
    # Manipulated variables
    m.beta.STATUS = 0  # calculated by the optimizer
    m.beta.FSTATUS = 0 # use measured value
    m.beta.DCOST = 0.0 # Delta cost penalty for MV movement
    m.beta.UPPER = 1.0 # Upper bound
    m.beta.LOWER = 0.0 # Lower bound
    m.beta.MV_STEP_HOR = 1
    m.beta.value = 0
    
    # Controlled variables
    m.Tb.STATUS = 1        # drive to set point
    m.Tb.FSTATUS = 0       # receive measurement
    m.Tb.SPHI = 17.5         # set point high level
    m.Tb.SPLO = 16.5         # set point low level
    m.Tb.WSPHI = 100       # set point high priority
    m.Tb.WSPLO = 100       # set point low priority
    
    T_MEAS = 20
    m.Tb.value = T_MEAS
    m.bias.value = T_MEAS - m.T.value[0]
    
    m.options.SOLVER = 3
    m.solve(disp=True)
    
    if m.options.APPSTATUS == 1:
       # Retrieve new values
       beta  = m.beta.NEWVAL
    else:
       # Solution failed
       beta  = 0.0
    
    # Plot the results
    plt.figure(figsize=(8,3.5))
    plt.subplot(3,1,1)
    plt.plot(m.time,m.Tb.value,'r-',label=r'$T_{int}$')
    plt.legend(loc=1); plt.grid()
    plt.ylabel('Tin (°C)')
    plt.subplot(3,1,2)
    plt.plot(m.time,m.d.value,'g:',label=r'$T_{ext}$')
    plt.ylabel('Tex (°C)')
    plt.legend(loc=1); plt.grid()
    plt.subplot(3,1,3)
    plt.step(m.time,m.beta.value,'b--',label=r'$\beta$')
    plt.ylabel('optimal control')
    plt.xlabel('Time (sec)')
    plt.legend(loc=1); plt.grid()
    plt.tight_layout()
    plt.savefig('results_step.png',dpi=300)
    plt.show()
    

    Here is the final controller with a sample sinusoidal profile for the external temperature.

    controller action

    # Import library
    import numpy as np
    import pandas as pd
    import time
    from gekko import GEKKO
    from numpy import array
    import matplotlib.pyplot as plt
    
    total_minutes_in_day = 24 * 60
    interval = 5  # minutes
    num_points = total_minutes_in_day // interval
    time_points = np.linspace(0, total_minutes_in_day, num_points)
    # Générer les points de temps
    T_externel = 8.5 + 11.5 * np.sin((2 * np.pi / total_minutes_in_day) * time_points - (np.pi / 2))
    
    #%% parameters
    p = {'a': array([[ 0.65059382],
           [-0.03027903],
           [ 0.06080624],
           [-0.02233488],
           [ 0.06877417],
           [ 0.37139369],
           [-0.22535397],
           [-0.01737772],
           [-0.03430238],
           [ 0.02420281]]),
         'b': array([[[ 3.71568256e-01, -8.09407185e+00],
            [-7.62375332e-02,  4.97505860e+00],
            [-1.34657379e-01, -1.56110911e+00],
            [ 7.32255512e-02,  1.89452386e+00],
            [-2.50583697e-02, -1.68565102e+00],
            [-2.06554465e-01,  2.03970469e+00],
            [ 1.50625153e-01,  8.79190526e-01],
            [-9.20974935e-03, -1.34505749e-01],
            [ 1.58409954e-03, -6.32745151e-01],
            [ 4.79076587e-03,  1.21199721e-01]]]),
         'c': array([2.08698176])}
    
    m = GEKKO(remote=False)
    m.y = m.Array(m.CV,1)
    m.u = m.Array(m.MV,2)
    m.arx(p,m.y,m.u)
    
    # rename CVs
    m.T = m.y[0]
    
    # rename MVs
    m.beta = m.u[1]
    
    # distrubance
    m.d = m.u[0]
    
    # distrubance and parametres
    m.d.value = 0
    
    m.bias = m.Param(0)
    m.Tb = m.CV()
    m.Equation(m.Tb==m.T+m.bias)
    
    # steady state initialization
    m.options.IMODE = 1
    m.solve(disp=False)
    
    m.d.value = T_externel
    
    m.options.IMODE   = 6 # MPC
    m.options.CV_TYPE = 1 # the objective is an l1-norm (region)
    m.options.NODES   = 2 # Collocation nodes
    m.options.SOLVER  = 3 # IPOPT
    m.time = time_points # step time = 300s
    
    # Manipulated variables
    m.beta.STATUS = 1  # calculated by the optimizer
    m.beta.FSTATUS = 0 # use measured value
    m.beta.DCOST = 1e-4 # Delta cost penalty for MV movement
    m.beta.UPPER = 1.0 # Upper bound
    m.beta.LOWER = 0.0 # Lower bound
    m.beta.MV_STEP_HOR = 1
    m.beta.value = 0
    
    # Controlled variables
    m.Tb.STATUS = 1        # drive to set point
    m.Tb.FSTATUS = 0       # receive measurement
    m.Tb.SPHI = 19         # set point high level
    m.Tb.SPLO = 17         # set point low level
    m.Tb.WSPHI = 100       # set point high priority
    m.Tb.WSPLO = 100       # set point low priority
    
    T_MEAS = 20
    m.Tb.value = T_MEAS
    m.bias.value = T_MEAS - m.T.value[0]
    
    m.options.SOLVER = 3
    m.solve(disp=True)
    
    if m.options.APPSTATUS == 1:
       # Retrieve new values
       beta  = m.beta.NEWVAL
    else:
       # Solution failed
       beta  = 0.0
    
    # Plot the results
    plt.figure(figsize=(8,3.5))
    plt.subplot(3,1,1)
    plt.plot(m.time,m.Tb.value,'r-',label=r'$T_{int}$')
    plt.plot([0,m.time[-1]],[m.Tb.SPHI,m.Tb.SPHI],'k--',label='Upper Bound')
    plt.plot([0,m.time[-1]],[m.Tb.SPLO,m.Tb.SPLO],'k--',label='Lower Bound')
    plt.legend(loc=1); plt.grid()
    plt.ylabel('Tin (°C)')
    plt.subplot(3,1,2)
    plt.plot(m.time,m.d.value,'g:',label=r'$T_{ext}$')
    plt.ylabel('Tex (°C)')
    plt.legend(loc=1); plt.grid()
    plt.subplot(3,1,3)
    plt.step(m.time,m.beta.value,'b--',label=r'$\beta$')
    plt.ylabel('optimal control')
    plt.xlabel('Time (sec)')
    plt.legend(loc=1); plt.grid()
    plt.tight_layout()
    plt.savefig('results_mpc.png',dpi=300)
    plt.show()