pythonoptimizationnonlinear-optimizationgekko

Hysteresis modelling as a control constraint for MPC in Python GEKKO

I am trying to introduce a hysteresis constraint in an MPC optimization problem for control signal dispatch using Python GEKKO. This has become a daunting task as I am unable to transform the following problem into equations that GEKKO understands. The problem:

An example of what I'm trying to do

Where:

engine = GEKKO(remote = False)

control = engine.Param(value = control_signal)

key        = engine.MV(value = 0)
key.STATUS =  1
key.LOWER  = -1
key.UPPER  =  1

engine.Equation(hysteresis_equation(key))

Manipulated variables in this case are dispatch percentiles of the control signal called 'key' which will influence the problem dynamics.

Where the key is the manipulated variable and hysteresis_equation is a function of the key value that should emulate a time dependent hysteresis. I have not given more details because there is no point, the problem resides in the implmentation of a non-linear hysteresis constraint in a GEKKO model.

I have tried looking at binary variables, however, I do not understand how to get them to change value throughout the optimization using GEKKO.

Tried calling an external function that returns True or False is not supported and yields @error: Equation Definition Equation without an equality (=) or inequality (>,<)

I have also tried introducing booleans in an equation that resembles power == (can_on * key + (1-can_on) *b0 + can_off * 0 + (1-can_off) * key)/2. The booleans are controlled variables in the hysteresis_equation and are set to 1 or 0 depending on the state of hysteresis, but are not GEKKO variables.

Thank you in advance for your help.

Solution

  • The easiest way to prevent the controller from turning off or on MVs too frequently is to use the u.MV_STEP_HOR option to specify the number of steps that it must be constant before it can move again. More details are in the documentation.

    Implementing a time-based condition is more complicated, but still possible. Below is an example script that has a minimum on-time of 2 cycles and a minimum off-time of 3 cycles. It uses logical conditions that switch the on_allow and off_allow states that control the inequality constraints on the MV derivative.

    On-time and off-time constraints

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

m = GEKKO()
m.time = np.linspace(0,10,11)

min_on_time = 2
min_off_time = 3

# Manipulated variable
u = m.MV(0, lb=0, ub=1, integer=True)
uc = m.Var(0); m.Equation(uc==u) # copy variable
u.STATUS = 1  # allow optimizer to change

is_on = m.Intermediate(u) # m.if3(u-0.5,0,1)
is_off = m.Intermediate(1-is_on)

on_time,off_time = m.Array(m.Var,2,value=0)
m.Equation(on_time.dt()==is_on)
m.Equation(off_time.dt()==is_off)

on_td,off_td = m.Array(m.Var,2,value=0)
m.delay(on_time,on_td,min_on_time+1)
m.delay(off_time,off_td,min_off_time+1)

off_allow = m.if3(on_time -on_td -min_on_time +0.5,0,1)
on_allow  = m.if3(off_time-off_td-min_off_time+0.5,0,1)

m.Equation(off_allow *(uc.dt()+1)+ (1-off_allow) * uc.dt()>=0)
m.Equation( on_allow *(uc.dt()-1)+ (1-on_allow)  * uc.dt()<=0)

# Process model
x = m.Var(value=0,lb=0,ub=5)
m.Equation(5*x.dt() == -x + 3*u)

# Objective
m.Minimize((x-2)**2)

m.options.IMODE = 6 # control
m.options.SOLVER = 1
m.solve(disp=True)

plt.figure(figsize=(6,3.5))
plt.subplot(3,1,1)
plt.step(m.time,u.value,'b-',label='u Profile')
plt.plot(m.time,x.value,'r--',label='x Response')
plt.legend(); plt.grid()
ax = plt.gca()
ax.xaxis.set_major_locator(MaxNLocator(integer=True))

plt.subplot(3,1,2)
plt.plot(m.time,on_allow.value,'r--',label='on allow')
plt.plot(m.time,off_allow.value,'k:',label='off allow')
ax = plt.gca()
ax.xaxis.set_major_locator(MaxNLocator(integer=True))
plt.legend(); plt.grid()

plt.subplot(3,1,3)
plt.plot(m.time,on_time.value,'r--',label='on time')
plt.plot(m.time,off_time.value,'k:',label='off time')
ax = plt.gca()
ax.xaxis.set_major_locator(MaxNLocator(integer=True))
plt.legend(); plt.grid()

plt.xlabel('Time'); plt.tight_layout()
plt.savefig('results.png',dpi=300)
plt.show()