Search code examples
gekko

How to add MV overshoot as a constraint in optimal PID parameters estimation problem?

This question is related to this my previous question how to formulate the problem of finding the optimal PID paramters in gekko?

I have successfully estimated the optimal PID parameters by minimizing the IAE based on the answer to the above question. However, I would like to add a constraint of MV overshoot in the objective function. Below is the pseudo code to add the constraint that the MV overshoot should be less than 10 % .

MV_Overshoot=  Max_MV - SteadyState_MV 
MV_Overshoot < 0.1*Steady_State MV  # constraint

However, I am stuck with the 2 questions below

  1. How to find the steady state_MV ?
  2. How to find the maximum MV from MV array ?

or is there any better way to calculate the MV overshoot ?

Solution

  • Solve a first step with a steady state calculation to determine the OP (MV) and PV values.

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

# Process model
Kp = 0.5                    # process gain
tauP = 10.0                 # process time constant

# pre-solve for steady state values
m = GEKKO()
OP = m.Var(value=0.0)       # controller output (MV)
PV = m.Var(value=0.0)       # process variable
SP = m.Param(value=5.0)     # final set point
m.Equation(SP==PV)          # assume PV reaches set point
m.Equation(tauP*PV.dt() + PV == Kp*OP)
m.options.IMODE = 1
m.solve(disp=False)

print('Steady state target values')
print('OP: ' + str(OP.value[0]))
print('PV: ' + str(PV.value[0]))
print('SP: ' + str(SP.value[0]))
PV_ss = PV.value[0]
OP_ss = OP.value[0]

    Solve a second step to optimize the tuning parameters to stay under the requested overshoot targets. Overshoot is typically defined in terms of the distance above the setpoint, not the amount over the MV target. Both are included in the example.

    Optimize

    # dynamic optimization of controller tuning
m = GEKKO()
tf = 40
m.time = np.linspace(0,tf,2*tf+1)
step = np.zeros(2*tf+1)
step[3:40] = 0.0
step[40:]  = 5.0

# Controller model
Kc = m.FV(15.0,lb=0,ub=100)   # controller gain
tauI = m.FV(100.0,lb=0.01,ub=100) # controller reset time
tauD = m.FV(0.0,lb=0,ub=100)  # derivative constant
Kc.STATUS = 1; tauI.STATUS = 1; tauD.STATUS = 1

OP_0 = m.Const(value=0.0)   # OP bias
OP = m.Var(value=0.0,ub=OP_ss*1.1) # controller output
                            #  with overshoot bound on OP (MV)
PV = m.Var(value=0.0,ub=PV_ss*1.1) # process variable
                            #  with overshoot bound on PV
SP = m.Param(value=step)    # set point
Intgl = m.Var(value=0.0)    # integral of the error
err = m.Intermediate(SP-PV) # set point error
m.Equation(Intgl.dt()==err) # integral of the error
m.Equation(OP == OP_0 + Kc*err + (Kc/tauI)*Intgl - PV.dt())
m.Equation(tauP*PV.dt() + PV == Kp*OP)
m.Minimize((SP-PV)**2)
m.options.IMODE=6
m.solve(disp=False)

print('Tuning values')
print('Kc: ' + str(Kc.value[0]))
print('tauI: ' + str(tauI.value[0]))
print('tauD: ' + str(tauD.value[0]))

plt.figure()
plt.subplot(2,1,1)
plt.plot([0,max(m.time)],[1.1*OP_ss,1.1*OP_ss],\
         'b-',label='Upper bound')
plt.plot(m.time,OP.value,'b:',label='OP')
plt.ylabel('Output')
plt.legend()
plt.subplot(2,1,2)
plt.plot(m.time,SP.value,'k-',label='SP')
plt.plot(m.time,PV.value,'r--',label='PV')
plt.plot([0,max(m.time)],[1.1*PV_ss,1.1*PV_ss],\
         'b-',label='Upper bound')
plt.xlabel('Time (sec)')
plt.ylabel('Process')
plt.legend()
plt.show()