Can Gekko change the value of a parameter multiple times during a Single simulation?

Question

I need to numerically solve the following equation set, which is a circuit model. I am using Gekko in Python to do this. Here is my equation: DAE

I need to do Simultaneous Simulation, so I am using IMODE=4, and here is my condensed code:

# Circuit parameters define
# ---

# Gekko setting
m = GEKKO(remote=True)
m.options.IMODE = 4  
m.options.COLDSTART = (
    2  # Cold start model: 0=warm start, 1=cold start, 2=decompose problem
)
m.options.REDUCE = 3  # Number of pre-processing cycles to identify equations or variables to eliminate

# Initial Values
# phin_dot means first-order derivative of phin
# phin_ddot means second-order derivative of phin
phi1 = m.Var(value=0)
phi1_dot = m.Var(value=0)  
phi1_ddot = m.Var(value=0)  
phi2 = m.Var(value=0)
phi2_dot = m.Var(value=0)
phi2_ddot = m.Var(value=0)
phi4 = m.Var(value=0)
phi4_dot = m.Var(value=0)
phi4_ddot = m.Var(value=0)
phi5 = m.Var(value=0)
phi5_dot = m.Var(value=0)
phi5_ddot = m.Var(value=0)
phi6 = m.Var(value=0)
phi6_dot = m.Var(value=0)
phi6_ddot = m.Var(value=0)
phi7 = m.Var(value=0)
phi7_dot = m.Var(value=0)
phi7_ddot = m.Var(value=0)
phi9 = m.Var(value=0)
phi9_dot = m.Var(value=0)

# Define DAE
m.Equation(
    [
        # Get derivatives
        phi1_dot == phi1.dt(),
        phi1_ddot == phi1_dot.dt(),
        phi2_dot == phi2.dt(),
        phi2_ddot == phi2_dot.dt(),
        phi4_dot == phi4.dt(),
        phi4_ddot == phi4_dot.dt(),
        phi5_dot == phi5.dt(),
        phi5_ddot == phi5_dot.dt(),
        phi6_dot == phi6.dt(),
        phi6_ddot == phi6_dot.dt(),
        phi7_dot == phi7.dt(),
        phi7_ddot == phi7_dot.dt(),
        phi9_dot == phi9.dt(),
        # DAE
        phi1_ddot
        == 1
        / C_J1
        * ((phi4 - phi1 - Phia1) / L1 - Ic1 * m.sin(phi1) - phi1_dot / R_J1),
        phi2_ddot
        == 1
        / C_J2
        * ((phi4 - phi2) / L2 - Ic2 * m.sin(phi2) - phi2_dot / R_J2),
        phi4_ddot
        == 1
        / Cd1
        * (
                (phi5 - phi4) / Ld1
                - (phi4 - phi1 - Phia1) / L1
                - (phi4 - phi2) / L2
        ),
        0
        == (phi9 - phi6 - Phia2) / L3
        + (phi9 - phi7) / L4
        - (phi5 - phi4) / Ld1,
        phi6_ddot - phi5_ddot
        == 1
        / C_J3
        * (
                (phi9 - phi6 - Phia2) / L3
                - Ic3 * m.sin(phi6 - phi5)
                - (phi6_dot - phi5_dot) / R_J3
        ),
        phi7_ddot - phi5_ddot
        == 1
        / C_J4
        * (
                (phi9 - phi7) / L4
                - Ic4 * m.sin(phi7 - phi5)
                - (phi7_dot - phi5_dot) / R_J4
        ),
        0 == (phi9 - phi6 - Phia2) / L3 + (phi9 - phi7) / L4 - ib,
    ]
)

# Solve equation
m.solve(disp=False)

# Plot results
# ---

Note that ib is a control parameter of this equation, not an unknown. The question is that I need to change ib many times during a Single simulation, in other words, ib is a vector [ib1, ib2, ..., ibn], and every ib in this vector will stay the same for the same amount of time, as shown below. I don't know how to deal with ib with Gekko, I have looked for the Gekko's documentation and Google, but still don't know how to implement this. Could someone help me, can Gekko do this?

ib will change with time

Answer 1

Gekko tutorials demonstrate how to give inputs to the model that are time-varying. Create a new parameter u with a NumPy array or Python list that has the changing values ux.

u = m.Param(value=ux)

A modification of problem 4 has multiple steps.

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

m = GEKKO()    # create GEKKO model
m.time = np.linspace(0,40,401) # time points

# create GEKKO parameter (steps)
ux = np.zeros(401)
for i,ui in enumerate(ux):
    if i==0:
        ux[i]=0
    else:
        if i%50==0:
            ux[i] = ux[i-1]+1
        else:
            ux[i] = ux[i-1]
ux[-100:] = ux[-101] # hold the last 100 constant
u = m.Param(value=ux)

x = m.Var(0.0) 
m.Equation(2*x.dt()==-x+0.5*u) 

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

plt.plot(m.time,u,'g:',label='u(t)')
plt.plot(m.time,x,'b-',label='x(t)')
plt.ylabel('values'); plt.xlabel('time'); plt.legend()
plt.show()

For long time horizons or a large model, try using m.options.IMODE=7 for sequential simulation to speed-up the code. Set disp=False to avoid printing the large output.