Dynamic Control Tuning

Dynamic controller tuning is the process of adjusting certain objective function terms to give more desirable solutions. As an example, a dynamic control application may either exhibit too aggressive manipulated variable movement or be too sluggish during set-point changes. Tuning is the iterative process of finding acceptable values that work over a wide range of operating conditions.

Common Tuning Parameters for MPC

Tuning typically involves adjustment of objective function terms or constraints that limit the rate of change (DMAX), penalize the rate of change (DCOST), or set absolute bounds (LOWER and UPPER). Measurement availability is indicated by the parameter (FSTATUS). The optimizer can also include (1=on) or exclude (0=off) a certain manipulated variable (MV) or controlled variable with STATUS. Below are common application, MV, and CV tuning constants that are adjusted to achieve desired model predictive control performance.

  • Application tuning
    • IMODE = 6 or 9 for model predictive control
    • CV_TYPE = CV type with 1=l_1-norm, 2=squared error
    • DIAGLEVEL = diagnostic level (0-10) for solution information
    • MAX_ITER = maximum iterations
    • MAX_TIME = maximum time before stopping
    • MV_TYPE = Set default MV type with 0=zero-order hold, 1=linear interpolation
    • SOLVER
      • 0=Try all available solvers
      • 1=APOPT (MINLP, Active Set SQP)
      • 2=BPOPT (NLP, Interior Point, SQP)
      • 3=IPOPT (NLP, Interior Point, SQP)
  • Manipulated Variable (MV) tuning
    • COST = (+) minimize MV, (-) maximize MV
    • DCOST = penalty for MV movement
    • DMAX = maximum that MV can move each cycle
    • FSTATUS = feedback status with 1=measured, 0=off
    • LOWER = lower MV bound
    • MV_TYPE = MV type with 0=zero-order hold, 1=linear interpolation
    • STATUS = turn on (1) or off (0) MV
    • UPPER = upper MV bound
  • Controlled Variable (CV) tuning
    • COST = (+) minimize CV, (-) maximize CV
    • FSTATUS = feedback status with 1=measured, 0=off
    • SP = set point with CV_TYPE = 2
    • SPLO = lower set point with CV_TYPE = 1
    • SPHI = upper set point with CV_TYPE = 1
    • STATUS = turn on (1) or off (0) CV
    • TAU = reference trajectory time-constant
    • TR_INIT = trajectory type, 0=dead-band, 1,2=trajectory
    • TR_OPEN = opening at initial point of trajectory compared to end

There are several ways to change the tuning values. Tuning values can either be specified before an application is initialized or while an application is running. To change a tuning value before the application is loaded, use the apm_option() function such as the following example to change the lower bound in MATLAB or Python for the MV named u.

 apm_option(server,app,'u.LOWER',0)

The upper and lower deadband targets for a CV named y are set to values around a desired set point of 10.0. In this case, an acceptable range for this CV is between 9.5 and 10.5.

 sp = 10.0
 apm_option(server,app,'y.SPHI',sp+0.5)
 apm_option(server,app,'y.SPHI',sp-0.5)

Application constants are modified by indicating that the constant belongs to the group nlc. IMODE is adjusted to either solve the MPC problem with a simultaneous (6) or sequential (9) method. In the case below, the application IMODE is changed to simultaneous mode.

 apm_option(server,app,'nlc.IMODE',6)

Exercise 1

Objective: Design a model predictive controller with one manipulated variable and two controlled variables with competing objectives that cannot be simultaneously satisfied. Tune the controller to achieve best performance. Estimated time: 2 hours.

Use the following system of linear differential equations for this exercise by placing the model definition in the file myModel.apm (APMonitor) or in a Python script (Gekko).

APMonitor Model

 Parameters
  u
 Variables
  x
  y
 Equations
  5 * $x = -x + 1.5 * u
  15 * $y = -y + 3.0 * u

Gekko Model

import numpy as np
from gekko import GEKKO
m = GEKKO()
t = np.linspace(0,15,76)
m.time = t # 0.2 cycle time
u = m.MV(1,name='u')
x = m.CV(name='x')
y = m.CV(name='y')
m.Equations([5*x.dt() ==-x+1.5*u,\
             15*y.dt()==-y+3.0*u])

In this case, the parameter u is the manipulated variable and x and y are the controlled variables. It is desired to maximize x and maintain values between 9 and 10. It is desired to maintain values of y between 2 and 7. For the first 10 minutes, the priority is to maintain the range for y and following this time period, it is desired to track the range for x. Tune the controller to meet these objectives while minimizing MV movement.

Solution

import numpy as np
from gekko import GEKKO
import matplotlib.pyplot as plt
import json
from IPython import display
from IPython import get_ipython

# detect if IPython is running
ipython = get_ipython() is not None

m = GEKKO(remote=False)
t = np.linspace(0,15,31)
m.time = t # 0.5 cycle time
u = m.MV(1,name='u')
x = m.CV(name='x')
y = m.CV(name='y')
m.Equations([5*x.dt() ==-x+1.5*u,\
             15*y.dt()==-y+3.0*u])
m.options.IMODE = 6
m.options.NODES = 3

# MV tuning
u.DCOST   = 0.1
u.DMAX    = 1.5
u.STATUS  = 1
u.FSTATUS = 0
u.UPPER   = 10
u.LOWER   = 0

# CV tuning
x.STATUS  = 1
x.FSTATUS = 0
x.SPHI    = 10
x.SPLO    = 9
x.TAU     = 2
x.TR_INIT = 1
x.TR_OPEN = 5
x.WSPHI   = 100
x.WSPLO   = 100

y.STATUS  = 1
y.FSTATUS = 0
y.SPHI    = 7
y.SPLO    = 2
y.TAU     = 0
y.TR_INIT = 0
y.TR_OPEN = 1
y.WSPHI   = 200
y.WSPLO   = 50

# Make an MP4 animation?
make_mp4 = True
if make_mp4:
    import imageio  # required to make animation
    import os
    try:
        os.mkdir('./figures')
    except:
        pass

# plot solution
plt.figure(figsize=(8,5))
plt.ion()
plt.show()

for i in range(40):
    m.solve(disp=False)

    # update time
    tm = m.time+i*0.5

    # x is more important at t=10
    if tm[0]>=10:
        x.WSPLO = 500
        x.WSPHI = 500

    # get trajectory information
    with open(m.path+'//results.json') as f:
        z = json.load(f)

    if ipython:
        display.clear_output(wait=True)
    else:
        plt.clf()

    plt.figure(1,figsize=(8,5))
    plt.subplot(3,1,1)
    plt.plot(tm,u.value,'k-',lw=2,label='u')
    plt.ylim([-1,11])
    plt.legend(loc=1)
    plt.ylabel('MV')

    plt.subplot(3,1,2)
    plt.plot(tm,x.value,'b-',lw=2,label='x')
    plt.plot(tm,z['x.tr_hi'],'k--',label='x traj')
    plt.plot(tm,z['x.tr_lo'],'k--')
    plt.legend(loc=1)

    plt.subplot(3,1,3)
    plt.plot(tm,y.value,'r-',lw=2,label='y')
    plt.plot(tm,z['y.tr_hi'],'k--',label='y traj')
    plt.plot(tm,z['y.tr_lo'],'k--')
    plt.legend(loc=1)
    plt.draw()
    if make_mp4:
        filename='./figures/plot_'+str(i+10000)+'.png'
        plt.savefig(filename)
    plt.pause(0.01)

# generate mp4 from png figures
if make_mp4:
    images = []
    for i in range(40):
        filename='./figures/plot_'+str(i+10000)+'.png'
        images.append(imageio.imread(filename))
    try:
        imageio.mimsave('results.mp4', images)
    except:
        print('Error: Install mp4 generator with "pip install imageio-ffmpeg"')

Exercise 2 — Surge Tank Level and Density Control

Objective: Design and tune an MPC that regulates (1) tank level within safe bounds and (2) product density to a target range despite feed-density disturbances. Use dilution-water flow as the single MV. Estimated time: 2–3 hours.

Process Description A cylindrical surge tank blends a dense slurry with dilution water. The outflow goes to a downstream unit that requires stable level and density.

  • Geometry
    • Tank height: H = 3.0 m
    • Tank diameter: D = 2.0 m
    • Cross-sectional area: A = πD²/4
  • Flows (m³/min)
    • Q_W: dilution-water inflow (MV)
    • Q_S: slurry inflow (disturbance)
    • Q_P: product outflow
  • States
    • Mass M (t) and volume V (m³)
    • Level (%): L = 100·V/(A·H)
    • Product density (kg/L): ρ_P = M/V
    • Water density ρ_W = 1.0 kg/L, slurry density ρ_S ≈ 2.0 kg/L + disturbances

Nonlinear Plant Model (mass/volume balances)

 dM/dt = Q_W·ρ_W + Q_S·ρ_S − Q_P·ρ_P
 dV/dt = Q_W + Q_S − Q_P
 L = 100·V/(A·H),  ρ_P = M/V

When L ≤ 0% or L ≥ 100%, hold dV/dt = 0 to avoid unphysical overflow/emptying during simulation.

Internal MPC Model A simplified internal model suitable for control design (you may keep the nonlinear balances as the “plant” simulator):

  • Level integrator (relative to Q_W* where * is steady state): (A·H)·(dL/dt) = 100·(Q_W − Q_W*)
  • First-order density response to dilution water: τ_model·(dρ_P/dt) = −(ρ_P − ρ_P*) + K_p·(Q_W − Q_W*) with K_p ≈ −0.24 (kg·L⁻¹)/(m³·min⁻¹), τ_model ≈ 0.8 min

Controller Design

  • MV: dilution-water flow Q_W
    • Bounds: 0.0 ≤ Q_W ≤ 1.6 m³/min (adjust as needed)
    • Movement limits: DMAX ≈ 0.1–0.3 m³/min per move
    • Movement penalty: DCOST ≈ 0.02–0.10
  • CVs:
    • Level L (%): Safety-critical. Keep within 20–80% at all times.
      • Use deadband tracking (CV_TYPE=1) with SPLO=20, SPHI=80.
      • Heavier penalties on violating bounds (e.g., WSPLO=WSPHI ≫ density weights).
    • Product density ρ_P (kg/L): Quality target 1.60 ± 0.05.
      • Deadband tracking with SPLO=1.55, SPHI=1.65 (optionally a small buffer inside these hard limits).
      • Moderate weight relative to level.
  • Priority / Importance
    • Level (safety) > Density (quality) > MV movement.
    • Example weights: L.WSPLO=L.WSPHI=10 (or higher if needed); ρ_P.WSPLO=ρ_P.WSPHI=1.
  • Timing
    • Controller interval: 1.0 min; prediction horizon ≈ 10–15 min; NODES≈2–3.

Disturbance Scenarios

  1. Random-walk changes in ρ_S with σ ≈ 0.018 kg/L per step (every 0.2 min).
  2. Optional steps: at t = 30 min, ρ_S → 2.10; at t = 60 min, ρ_S → 1.90.

Tasks

  1. Implement the internal MPC model (level integrator + first-order density).
  2. Tune DMAX, DCOST, and CV weights to:
    1. keep L within 20–80% at all times (no overflow/empty),
    2. keep ρ_P within 1.60 ± 0.05 after disturbances (settle < ~3–5 min),
    3. avoid excessive Q_W movement.

# --- Imports ---
import math as m
import numpy as np
from scipy.optimize import fsolve
from scipy.integrate import odeint
from gekko import GEKKO
from matplotlib import pyplot as plt
import random

# --- Geometry ---
tank_height  = 3.0  # m
tank_diameter= 2.0  # m
tank_area    = (m.pi/4.0) * (tank_diameter**2)  # m^2

# --- Nonlinear plant model (for simulation) ---
def plant(states, t, Q_W, Q_S, rho_S, Q_P):
    M, V = states
    L = 100.0 * V / (tank_height * tank_area)
    rho_P = M / V
    rho_W = 1.0  # kg/L
    dMdt = (Q_W * rho_W + Q_S * rho_S - Q_P * rho_P)
    dVdt = Q_W + Q_S - Q_P
    # prevent unphysical empty/overflow dynamics
    if L <= 0.0 or L >= 100.0:
        dVdt = 0.0
    return [dMdt, dVdt]

# --- Nominal targets ---
ss_rho_S = 2.0   # kg/L (slurry)
ss_Q_S   = 1.5   # m^3/min (slurry inflow)
ss_rho_P = 1.60  # kg/L (product target)
ss_level = 50.0  # % level

# --- Solve steady state (Q_W*, Q_P*, M*, V*) ---
def ss_residuals(z):
    Q_W, Q_P, M, V = z
    dMdt, dVdt = plant([M, V], 0.0, Q_W, ss_Q_S, ss_rho_S, Q_P)
    eq1 = dMdt
    eq2 = dVdt
    eq3 = ss_rho_P - (M / V)                      # density target
    eq4 = ss_level - 100.0*V/(tank_height*tank_area)  # level target
    return [eq1, eq2, eq3, eq4]

ss_guess = [1.0, 1.0, 10.0, 10.0]
ss_Q_W, ss_Q_P, ss_M, ss_V = fsolve(ss_residuals, ss_guess)

# --- MPC model (internal) ---
m = GEKKO(remote=False)
m.options.IMODE = 6            # MPC
m.options.CV_TYPE = 1          # deadband tracking
m.options.DIAGLEVEL = 0

# time grid (controller executes every 1.0 min; dt = 0.2 min)
tmax_gekko = 15.0
dt_gekko   = 0.2
m.time = np.linspace(0, tmax_gekko, int(tmax_gekko/dt_gekko)+1)
m.options.MV_STEP_HOR = int(1.0/dt_gekko)  # 1-min move interval

# CVs
rho_P_cv = m.CV(name='rho_P'); rho_P_cv.STATUS=1; rho_P_cv.FSTATUS=1
rho_P_cv.WSPHI = 1; rho_P_cv.WSPLO = 1
density_db = 0.05
rho_P_cv.SPHI = ss_rho_P + density_db - 0.035  # inside hard limits
rho_P_cv.SPLO = ss_rho_P - density_db + 0.035
rho_P_cv.value = ss_rho_P

level_cv = m.CV(name='level'); level_cv.STATUS=1; level_cv.FSTATUS=1
level_cv.SPHI = 80.0; level_cv.SPLO = 20.0
level_cv.WSPHI = 10; level_cv.WSPLO = 10
level_cv.value = ss_level

# MV
QW = m.MV(name='QW'); QW.STATUS=1; QW.DCOST=0.05; QW.LOWER=0.0
QW.value = ss_Q_W

# Internal model parameters
Kp_model  = -0.24   # (kg/L)/(m^3/min)
tau_model = 0.8     # min

# Internal model equations (relative to steady state)
m.Equation((tank_area*tank_height)*level_cv.dt() == 100.0*(QW - ss_Q_W))
m.Equation(tau_model * rho_P_cv.dt() == -(rho_P_cv - ss_rho_P) + Kp_model*(QW - ss_Q_W))

# --- Figure scaffolding (optional) ---
Fig, axs = plt.subplots(4,1, figsize=(7, 7.2), sharex=True)
axs[0].axhline(level_cv.SPHI, linestyle='--'); axs[0].axhline(level_cv.SPLO, linestyle='--')
axs[0].set_ylim(0,100); axs[0].set_ylabel('Level (%)')

axs[1].axhline(ss_rho_P + density_db, linestyle='--')
axs[1].axhline(ss_rho_P - density_db, linestyle='--')
axs[1].set_ylim(1.45,1.75); axs[1].set_ylabel('Density (kg/L)')

axs[2].set_ylim(0.5,1.6); axs[2].set_ylabel('Q_W (m³/min)')
axs[3].set_ylim(1.75,2.3); axs[3].set_ylabel('ρ_S (kg/L)')
axs[-1].set_xlabel('Time (min)')

# --- Your tasks from here ---
# 1) Create a time loop, generate ρ_S disturbances (random walk),
# 2) at each control interval: set MEAS (level_cv.MEAS, rho_P_cv.MEAS), m.solve(disp=False),
# 3) integrate the plant by one step with the chosen QW.NEWVAL,
# 4) log/plot L, ρ_P, Q_W, ρ_S and verify specs.

```

Full Solution

# import required libraries
import math as m
import numpy as np
from scipy.optimize import fsolve
from gekko import GEKKO
from matplotlib import pyplot as plt
from scipy.integrate import odeint
import random

# define the surge tank dimensions

tank_height = 3.0 # m
tank_diameter = 2.0 # m

# Cross-sectional area of cylindrical tank (m²)
tank_area = (m.pi / 4) * (tank_diameter ** 2)

# define the non-linear model
def model(states, t, Q_W, Q_S, rho_S, Q_P):

    '''
    Computes mass and volume rate changes in a surge tank.
        '''


    # unpack the states
    mass, vol = states # mass in t and vol in m3

    # calculate the level %
    L = 100 * vol / (tank_height * tank_area)

    # Calculate density of material in tank (t/m3 equivalent to kg/L)
    rho_P = mass / vol

    # Density of dilution water
    rho_W = 1.0 # kg/L

    # Mass balance: input from dilution water and slurry minus output from product stream
    dmdt = (Q_W * rho_W + Q_S * rho_S - Q_P * rho_P) # t/min

    # Volume balance: sum of inflows minus outflow
    dvdt = Q_W + Q_S - Q_P # m3/min

    # prevent change in volume if tank is empty or overflowing
    if L <=0 or L >=100:
        dvdt = 0

    return [dmdt, dvdt]


# set the desired steady-state conditions
ss_rho_S = 2.0 # density of feed slurry (kg/L)
ss_rho_P = 1.6 # density of product (kg/L)
ss_Q_S = 1.5 # flow rate of slurry (m3/min)
ss_level = 50 # steady-state level of the tank (%)

# define a function used to solve for the steady-state

def steady_state(unknowns):

    """
    Computes the residuals for the steady-state mass and volume balance of a surge tank system.
    """


    Q_W, Q_P, mass, vol = unknowns # unknowns to solve for

    eqs = [0] * len(unknowns) # list of residuals with placeholder 0s

    derivs = model([mass, vol], 0, Q_W, ss_Q_S, ss_rho_S, Q_P)

    eqs[0] = derivs[0] # Mass derivative (should be zero at steady-state)
    eqs[1] = derivs[1] # Volume derivative (should be zero at steady-state)
    eqs[2] = ss_rho_P - mass / vol # Difference between target and calculated product density
    eqs[3] = ss_level - 100 * vol / (tank_height * tank_area) # Difference between target and calculated tank level

    return eqs

# Initial guess for the steady-state values:
guess = [1, 1, 10, 10]

# Solve the system of equations defined in the steady_state function
ss_sol = fsolve(steady_state, guess)

# Extract the steady-state dilution water flow rate
ss_Q_W = ss_sol[0]

# Extract the steady-state product flow rate
ss_Q_P = ss_sol[1]

# Extract the steady-state mass of the tank contents
ss_mass = ss_sol[2]

# Extract the steady-state volume of the tank contents
ss_vol = ss_sol[3]

##  CONTROLLER DEVELOPMENT ##

# create and configure gekko model
m = GEKKO(remote=False)
m.options.IMODE = 6
m.options.CV_TYPE = 1
# m.open_folder()
m.options.DIAGLEVEL = 0

tmax_gekko = 15
dt_gekko = 0.2
m.time = np.linspace(0,tmax_gekko,int(tmax_gekko / dt_gekko) + 1) # gekko model time span

controller_interval = 1.0 # how often the controller compute new control actions? (min)

m.options.MV_STEP_HOR = int(controller_interval / dt_gekko)

# create density variable
rho_P_gekko = m.CV(name = 'rho_P')
rho_P_gekko.value = ss_rho_P
rho_P_gekko.STATUS = 1
rho_P_gekko.WSPHI = 1
rho_P_gekko.WSPLO = 1
rho_P_gekko.WMEAS = 0
rho_P_gekko.FSTATUS = 1
rho_P_gekko.COST = 0
# rho_P_gekko.SP = ss_rho_P

density_db = 0.05 # density deadband (kg/L)
density_Hi_limit = ss_rho_P + density_db
density_Lo_limit = ss_rho_P - density_db

density_db_buffer = 0.035 # distance between SPs and limits

rho_P_gekko.SPHI = density_Hi_limit - density_db_buffer
rho_P_gekko.SPLO = density_Lo_limit + density_db_buffer

# create level variable
level_gekko = m.CV(name = 'level')
level_gekko.value = ss_level
# level_gekko.SP = 50
level_gekko.SPHI = 80
level_gekko.SPLO = 20
level_gekko.STATUS = 1
level_gekko.WSPHI = 10 # heavy penalties for exceeding level limits
level_gekko.WSPLO = 10
level_gekko.FSTATUS = 1
level_gekko.COST = 0

# create dilution water variable
Q_W_gekko = m.MV(name = 'QW')
Q_W_gekko.value = ss_Q_W
Q_W_gekko.STATUS = 1
Q_W_gekko.DCOST = 0.05
Q_W_gekko.LOWER = 0.0
Q_W_gekko.COST = 0
Q_W_gekko.MV_STEP_HOR = 0 # use global step horizon

# define the 1st-order density model parameters
Kp_model = -0.24 # Steady-state gain in density for a change in water flow ((g/L) / (m3/min))
tau_model = 0.8 # approximated as +- residence time

# define the internal controller equations
m.Equation((tank_area * tank_height) * level_gekko.dt() == 100 * (Q_W_gekko - ss_Q_W))
m.Equation(tau_model * rho_P_gekko.dt() == -(rho_P_gekko - ss_rho_P) + Kp_model*(Q_W_gekko - ss_Q_W))


##  CREATE FIGURE FOR PLOTTING ##

Fig, axs = plt.subplots(4,1, figsize = (1*7, 4*1.8), sharex = True)

# Subplot 0: Tank level
axs[0].axhline(level_gekko.SPHI, color='red', linestyle='--')
axs[0].axhline(level_gekko.SPLO, color='red', linestyle='--')
axs[0].set_ylim(0, 100)
axs[0].set_ylabel('Tank level (%)')

level_line, = axs[0].plot([],[], color = 'blue')
level_line_pred, = axs[0].plot([],[], color = 'blue', linestyle = '--')
now_line0 = axs[0].axvline([0])

# Subplot 1: Product density
axs[1].axhline(density_Hi_limit, color='red', linestyle='--')
axs[1].axhline(density_Lo_limit, color='red', linestyle='--')
axs[1].set_ylim(1.45, 1.75)
axs[1].set_ylabel('Product \ndensity (kg.L$^{-1}$)')

density_line, = axs[1].plot([],[], color = 'blue')
density_line_pred, = axs[1].plot([],[], color = 'blue', linestyle = '--')
now_line1 = axs[1].axvline([0])

# Subplot 2: Dilution water flow rate
axs[2].set_ylim(0.5, 1.6)
axs[2].set_ylabel('Dilution water \nflow rate (m$^3$.min$^{-1}$)')

water_line, = axs[2].plot([],[], color = 'blue')
water_line_pred, = axs[2].plot([],[], color = 'blue', linestyle = '--')

now_line2 = axs[2].axvline([0])

# Subplot 3: Feed slurry density (FIXED)
axs[3].set_ylim(1.75, 2.3)
axs[3].set_ylabel('Feed slurry \n density (kg.L$^{-1}$)')

slurry_density_line, = axs[3].plot([],[], color = 'blue')
now_line3 = axs[3].axvline([0])

# Set common x-axis
axs[-1].set_xlabel('Time (min)')
axs[-1].set_xlim(0,90)

##  SIMULATION ##

# set the simulation time frame
tsim = 90 # min
dt_sim = dt_gekko
nt = int(tsim/dt_sim) + 1
tspan = np.linspace(0, tsim, nt)

# set the random slurry densities at each time step
random.seed(5)
increments = [random.gauss(0, 0.018) for _ in range(nt - 1)]
rho_S = [ss_rho_S] + [ss_rho_S + sum(increments[:i]) for i in range(1, nt)]

# create data storage lists with placeholder 0s
Q_W = [0] * nt

rho_P = [0] * nt
rho_P[0] = ss_rho_P

level = [0] * nt
level[0] = ss_level

# inital states from which to start the simulation
states = [ss_mass, ss_vol]

# initial time of the next controller calculation
next_solve_time = 0 # min

# simulation loop
for i in range(nt-1):

    t = tspan[i]

    if t >= next_solve_time:

        level_gekko.MEAS = level[i]
        rho_P_gekko.MEAS = rho_P[i]

        m.solve(disp = False)
        # input('')
        Q_calc = Q_W_gekko.NEWVAL

        next_solve_time += controller_interval

    # hold constant values for slurry and product flow rates
    Q_S = ss_Q_S
    Q_P = ss_Q_P

    # integrate the non-linear model by one time step to calculate the new states
    states = odeint(model, states, [t, t+dt_sim], args = (Q_calc, Q_S, rho_S[i], Q_P))[-1].tolist()

    mass = states[0]
    volume = states[1]

    # store the data for the time step
    Q_W[i] = Q_calc
    rho_P[i+1] = mass / volume
    level[i+1] = 100 * volume / (tank_height * tank_area)

    # plotting

    level_line.set_data(tspan[:i], level[:i])
    level_line_pred.set_data(np.linspace(next_solve_time, next_solve_time + tmax_gekko, len(level_gekko.PRED)-1), level_gekko.PRED[1:])
    density_line.set_data(tspan[:i], rho_P[:i])
    density_line_pred.set_data(np.linspace(next_solve_time, next_solve_time + tmax_gekko, len(rho_P_gekko.PRED)-1), rho_P_gekko.PRED[1:])
    water_line.set_data(tspan[:i], Q_W[:i])
    water_line_pred.set_data(np.linspace(next_solve_time, next_solve_time + tmax_gekko, len(Q_W_gekko.PRED)-1), Q_W_gekko.PRED[1:])
    slurry_density_line.set_data(tspan[:i], rho_S[:i])

    now_line0.set_xdata([t])
    now_line1.set_xdata([t])
    now_line2.set_xdata([t])
    now_line3.set_xdata([t])

    plt.pause(0.2)

Thanks to Vincent Roos for the Level and Density control exercise and Gekko solution.


Discussion Questions

  1. A hard upper bound on a controlled state and a soft constraint with a large violation penalty can produce the same closed-loop trajectory in nominal operation. Under what disturbance do they behave completely differently, and which failure mode do you prefer in a real plant?
  2. The squared-error objective (CV_TYPE=2) with a reference trajectory and the l1 deadband objective (CV_TYPE=1) both "track a setpoint." Which one keeps moving the MV when the CV sits inside a noise band, and what does that cost in actuator life?
  3. DMAX and DCOST both slow the manipulated variable. One is a constraint, one is a price. Which one can make the problem infeasible, which one degrades gracefully, and when would you use each?
  4. Increasing TAU on the reference trajectory and increasing DCOST on the MV both make the response gentler. How would you tell from closed-loop data which one a previous engineer used?

Generative AI Learning

Use these prompts with the tuning exercises (competing CVs, and the surge tank level/density problem). Direct the AI - the tuning tradeoffs are engineering judgment, which is the deliverable.

"Quiz me with 5 questions, one at a time, on MPC tuning in GEKKO/APMonitor: what DCOST and DMAX do to a manipulated variable and how they differ in failure mode, what TAU and TR_INIT define for a controlled variable's reference trajectory, how WSPHI/WSPLO with SPHI/SPLO implement an asymmetric deadband, what CV_TYPE selects and when the l1 form is preferred, and why one MV serving two CVs (level and density) cannot satisfy both when their required moves conflict. Grade my answers and list my misconceptions."
"Here is my surge tank MPC tuning: {MV: DCOST, DMAX; level CV: SPHI/SPLO, weights; density CV: setpoint/band, weights} and closed-loop plots described as {describe}. Attack it as a plant operations reviewer: find the disturbance scenario (feed density step, feed flow surge) where my weights let the tank overflow or run empty before quality is sacrificed, ask me to defend the priority ordering of level vs density, and make me state the largest feed-density disturbance my tuning survives - with the calculation, not a guess."

What to Turn In

Submit a short report (PDF, 1-2 pages) that curates your results into a demonstration of what you learned. You may use Generative AI to help write the report, but you must guide it to the correct visualizations, justifications, and assumptions. Answer these questions:

  1. Exercise 1: show the two-CV solution with the time-based priority switch, and explain which tuning values enforce the priority and how the MV movement was kept minimal.
  2. Exercise 2 (surge tank): include level and density plots against the feed disturbances, your final tuning table, and one sentence per parameter on why that value.
  3. Demonstrate one hard-vs-soft constraint comparison on the tank level: same disturbance, both formulations. What differs in the trajectories and in solver behavior?
  4. From the reviewer prompt: the disturbance scenario that broke (or nearly broke) your tuning, and what you changed in response.

Course Information

TCLab

Project

Applications

Exams

Modeling

Machine Learning

Estimation

Control

Reinforcement Learning

Related Courses

Admin