Proportional-Integral-Derivative (PID) Controller#

The IDAES framework contains a basic PID control implementation, which is described in this section.

Example#

The following code provides a full example for using the PIDController model. The flowsheet model is a tank with a valve on the inlet and outlet. The dynamics of the valves are assumed to be very fast and is not considered. The system has steam flowing through it and the inlet valve is used to maintain the tank pressure at a given setpoint. The final model solved demonstrates the controller with a step change in the inlet pressure. Results are saved to a file “pid_steam_tank_pressure.pdf.”

import pyomo.environ as pyo
from pyomo.network import Arc
from idaes.core import FlowsheetBlock, MaterialBalanceType
from idaes.models.unit_models import Heater, Valve
from idaes.models.properties import iapws95
from idaes.core.util.initialization import propagate_state
from idaes.models.control.controller import (
    PIDController,
    ControllerType,
    ControllerMVBoundType
)
import idaes.core.util.scaling as iscale
from idaes.core.solvers import get_solver
from idaes.core.util.plot import plot_grid_dynamic

def _valve_pressure_flow_cb(b):
    """Callback for the pressure flow relations for the valves"""
    umeta = b.control_volume.config.property_package.get_metadata().get_derived_units
    b.Cv = pyo.Var(
        initialize=0.1,
        doc="Valve flow coefficient",
        units=umeta("amount") / umeta("time") / umeta("pressure"),
    )
    b.Cv.fix()
    b.flow_var = pyo.Reference(b.control_volume.properties_in[:].flow_mol)
    b.pressure_flow_equation_scale = lambda x: x ** 2
    @b.Constraint(b.flowsheet().time)
    def pressure_flow_equation(b2, t):
        Po = b2.control_volume.properties_out[t].pressure
        Pi = b2.control_volume.properties_in[t].pressure
        F = b2.control_volume.properties_in[t].flow_mol
        Cv = b2.Cv
        fun = b2.valve_function[t]
        return F ** 2 == Cv ** 2 * (Pi ** 2 - Po ** 2) * fun ** 2

def _add_inlet_pressure_step(m, time=1, value=6.0e5):
    """Add an inlet pressure step change"""
    for t in m.fs.time:
        if t >= time:
            m.fs.valve_1.inlet.pressure[t].fix(value)

def create_model(time_set=None, time_units=pyo.units.s, nfe=5, tee=False):
    """Build and initialize the flowsheet model

           valve_1   +----+
  steam ----|><|-->--| ta |    valve_2
                     | nk |-----|><|--->--- steam
                     +----+
    """
    m = pyo.ConcreteModel(name="Dynamic Steam Tank with PID Control")
    # Add flowsheet
    m.fs = FlowsheetBlock(dynamic=True, time_set=time_set, time_units=time_units)
    # Add water property parameter block
    m.fs.prop_water = iapws95.Iapws95ParameterBlock(
        phase_presentation=iapws95.PhaseType.LG
    )
    # Add valve 1
    m.fs.valve_1 = Valve(
        dynamic=False,
        has_holdup=False,
        pressure_flow_callback=_valve_pressure_flow_cb,
        material_balance_type=MaterialBalanceType.componentTotal,
        property_package=m.fs.prop_water,
    )
    # Add heater model to represent a tank (close to bare control volume model).
    m.fs.tank = Heater(
        has_holdup=True,
        material_balance_type=MaterialBalanceType.componentTotal,
        property_package=m.fs.prop_water,
    )
    # Add valve 2
    m.fs.valve_2 = Valve(
        dynamic=False,
        has_holdup=False,
        pressure_flow_callback=_valve_pressure_flow_cb,
        material_balance_type=MaterialBalanceType.componentTotal,
        property_package=m.fs.prop_water,
    )
    # Add a controller
    m.fs.ctrl = PIDController(
        process_var=m.fs.tank.control_volume.properties_out[:].pressure,
        manipulated_var=m.fs.valve_1.valve_opening,
        calculate_initial_integral=True,
        mv_bound_type=ControllerMVBoundType.SMOOTH_BOUND,
        controller_type=ControllerType.PI,
    )
    # The control volume block doesn't assume equilibrium, so I'll make that
    # assumption here. I don't actually expect liquid to form but who knows?
    # The phase_fraction in the control volume is volumetric phase fraction.
    @m.fs.tank.Constraint(m.fs.time)
    def vol_frac_vap(b, t):
        return (
            b.control_volume.properties_out[t].phase_frac["Vap"]
            * b.control_volume.properties_out[t].dens_mol
            / b.control_volume.properties_out[t].dens_mol_phase["Vap"]
        ) == (b.control_volume.phase_fraction[t, "Vap"])

    # Connect the models
    m.fs.v1_to_tank = Arc(source=m.fs.valve_1.outlet, destination=m.fs.tank.inlet)
    m.fs.tank_to_v2 = Arc(source=m.fs.tank.outlet, destination=m.fs.valve_2.inlet)

    # Add the stream constraints
    pyo.TransformationFactory("network.expand_arcs").apply_to(m.fs)

    # Do DAE discretization
    pyo.TransformationFactory("dae.finite_difference").apply_to(
        m.fs,
        nfe=nfe,
        wrt=m.fs.time,
        scheme="BACKWARD"
    )

    # Fix the derivative variables to zero at time 0 (steady state assumption)
    m.fs.fix_initial_conditions()

    # Fix the input variables
    m.fs.valve_1.Cv.fix(0.001) # valve 1 flow coefficient
    m.fs.valve_2.Cv.fix(0.001) # valve 2 flow coefficient
    m.fs.ctrl.gain_p.fix(1e-6) # proportional gain
    m.fs.ctrl.gain_i.fix(1e-5) # integral error gain
    m.fs.ctrl.setpoint.fix(3e5) # setpoint is tank pressure of 300 kPa
    m.fs.ctrl.mv_ref.fix(0) # controller bias
    m.fs.ctrl.mv_lb = 0.0 # valve opening lower bound is 0
    m.fs.ctrl.mv_ub = 1.0 # value opening upper bound is 1
    m.fs.valve_1.inlet.enth_mol.fix(50000) # inlet fluid molar enthalpy
    m.fs.valve_1.inlet.pressure.fix(5e5) # inlet pressure of valve 1
    m.fs.valve_2.outlet.pressure.fix(101325) # valve 2 outlet pressure
    m.fs.valve_1.valve_opening.fix(1) # fix valve 1 opening to full open
    m.fs.valve_2.valve_opening.fix(1) # fix valve 2 opening to full open
    m.fs.tank.heat_duty.fix(0) # no heat transfer to tank
    m.fs.tank.control_volume.volume.fix(2.0) # 2 m^3 tank volume

    # Initialize the model
    solver = get_solver(options={"max_iter": 50})

    for t in m.fs.time:
        m.fs.valve_1.inlet.flow_mol[t] = 100  # initial guess on flow
    # simple initialize
    m.fs.valve_1.initialize()
    propagate_state(m.fs.v1_to_tank)
    m.fs.tank.initialize()
    propagate_state(m.fs.tank_to_v2)
    # Can't specify both flow and outlet pressure so free the outlet pressure
    # for initialization and refix it after.  Inlet flow gets fixed in init
    op = {}
    for t in m.fs.time:
        op[t] = pyo.value(m.fs.valve_2.outlet.pressure[t])
        m.fs.valve_2.outlet.pressure[t].unfix()
    m.fs.valve_2.initialize()
    for t in m.fs.time:
        m.fs.valve_2.outlet.pressure[t].fix(op[t])

    # Solve dynamic problem with deactivated controller
    m.fs.ctrl.deactivate()
    m.fs.valve_1.valve_opening.fix()
    solver.solve(m, tee=tee)

    # Solve dynamic problem with controller on
    m.fs.ctrl.activate()
    m.fs.valve_1.valve_opening.unfix()
    m.fs.valve_1.valve_opening[m.fs.time.first()].fix()
    solver.solve(m, tee=tee)

    # Return the model and solver
    return m, solver

# Create a model for the 0 to 12 sec time period
m, solver = create_model(time_set=[0, 12], nfe=30, tee=False)

# Add a step change in inlet pressure at t=6 from 500 kPa to 550 kPa
_add_inlet_pressure_step(m, time=6, value=5.5e5)

# Solve with step change
solver.solve(m, tee=False)

# Plot some results and save the figure to a file.
plot_grid_dynamic(
    x=m.fs.time,
    xlabel="time (s)",
    y=[
        m.fs.valve_1.valve_opening,
        m.fs.tank.control_volume.properties_out[:].pressure,
        m.fs.valve_1.control_volume.properties_in[:].pressure,
    ],
    ylabel=[
        "opening (fraction open)",
        "tank pressure (kPa)",
        "inlet pressure (kPa)",
    ],
    yunits=[
        None,
        pyo.units.kPa,
        pyo.units.kPa,
    ],
    cols=3,
    rows=1,
    to_file="pid_steam_tank_pressure.pdf"
)

assert abs(
  pyo.value(m.fs.valve_1.valve_opening[m.fs.time.last()]) - 0.61254) < 0.001

Class Documentation#

class idaes.models.control.controller.PIDController(*args, **kwds)#

PID controller model block. To use this the model must be dynamic.

Parameters:
  • rule (function) – A rule function or None. Default rule calls build().

  • concrete (bool) – If True, make this a toplevel model. Default - False.

  • ctype (class) –

    Pyomo ctype of the block. Default - pyomo.environ.Block

    Config args

    dynamic

    Indicates whether this model will be dynamic or not, default = useDefault. Valid values: { useDefault - get flag from parent (default = False), True - set as a dynamic model, False - set as a steady-state model.}

    has_holdup

    Indicates whether holdup terms should be constructed or not. Must be True if dynamic = True, default - False. Valid values: { useDefault - get flag from parent (default = False), True - construct holdup terms, False - do not construct holdup terms}

    process_var

    A Pyomo Var, Expression, or Reference for the measured process variable. Should be indexed by time. Slices are okay.

    manipulated_var

    A Pyomo Var, Reference to a Var, or something that can be used to construct a Reference to a Var for the controlled variable. The final Reference should be indexed by time. Slices are okay.

    mv_bound_type

    Type of bounds to apply to the manipulated variable output. If, bounds are applied, the model parameters mv_lb and mv_ub set the bounds. The default is ControllerMVBoundType.NONE. See the controller documentation for details on the mathematical formulation. The options are: ControllerMVBoundType.NONE no bounds, ControllerMVBoundType.SMOOTH_BOUND smoothed mv = min(max(mv_unbound, ub), lb), and ControllerMVBoundType.LOGISTIC logistic function to enforce bounds.

    calculate_initial_integral

    Calculate the initial integral term value if True

    controller_type

    Controller type. The default = ControllerType.PI and the options are: ControllerType.P Proportional, ControllerType.PI proportional and integral, ControllerType.PD proportional and derivative, and ControllerType.PID proportional, integral, and derivative

    antiwindup_type

    Type of antiwindup technique to use. Options are ControllerAntiwindupType.NONE, ControllerAntiwindupType.CONDITIONAL_INTEGRATION, and ControllerAntiwindupType.BACK_CALCULATION. See the controller documentation for details on the mathematical formulation.

    derivative_on_error

    Naive implementations of derivative action can cause large spikes in control when the setpoint is changed. One solution is to use the (negative) derivative of the process variable to calculate derivative action instead of using the derivative of setpoint error. If True, use the derivative of setpoint error to calculate derivative action. If False (default), use the (negative) derivative of the process variable instead.

  • initialize (dict) – ProcessBlockData config for individual elements. Keys are BlockData indexes and values are dictionaries with config arguments as keys.

  • idx_map (function) – Function to take the index of a BlockData element and return the index in the initialize dict from which to read arguments. This can be provided to override the default behavior of matching the BlockData index exactly to the index in initialize.

Returns:

(PIDController) New instance

class idaes.models.control.controller.PIDControllerData(component)[source]#

PID controller class.

build()[source]#

Build the PID block

Variables, Parameters, and Expressions#

Symbol

Name in Model

Description

\(y_{sp}(t)\)

setpoint[t]

Setpoint variable

\(y(t)\)

process_var[t]

Measured process variable (Reference)

\(u(t)\)

manipulated_var[t]

Manipulated variable (Reference)

\(K_p(t)\)

gain_p[t]

Controller gain (usually fixed)

\(K_i(t)\)

gain_i[t]

Integral gain (usually fixed)

\(K_d(t)\)

gain_d[t]

Derivative gain (usually fixed)

\(K_b(t)\)

gain_b[t]

Back-calculation gain (usually fixed)

\(u_{ref}\)

mv_ref[t]

Reference value of manipulated variable

\(e(t)\)

error[t]

Error expression or variable (setpoint - pv)

\(e_d(t) \quad\text{or}\; -\frac{dy(t)}{dt}\)

derivative_term[t]

Derivative of either error or negative pv

\(u_i(t)\)

mv_integral_component[t]

Portion of control from integral action

\(\frac{de_i(t)}{dt}\)

mv_integral_component_dot[t]

Derivative of integral term w.r.t. time

\(u_{ub}\)

mv_ub

Upper limit of output parameter

\(u_{lb}\)

mv_lb

Lower limit of output parameter

\(\epsilon\)

smooth_eps

Smooth min/max smoothing parameter

\(k_\text{log}\)

logistic_bound_k

Logistic bound steepness parameter

\(k_\text{cond}\)

conditional_integration_k

Conditional integration steepness parameter

Formulation#

Textbook PID#

The PID controller equations are given by the following equations (see the Variables, Parameters and Expressions section for a description of the variables).

The PID unbounded controller model is given by the following equation where the integral and derivative terms are optional.

\[u(t) = K_p e(t) + K_i e_i(t) + K_d e_d(t) + u_{ref}(t)\]

The error term is given by the following equation.

\[e(t) = y_{sp}(t) - y(t)\]

The integral error is defined as follows, where the initial condition of the integral error can be calculated from the manipulated variable initial condition to ensure that the initial conditions are consistent and avoid an initial jump in the manipulated variable value.

\[e_i(t) = e_i(t_0) + \int_{t_0}^t e(\tau) \text{d}\tau\]

The derivative of error is given below.

\[e_d(t) = \frac{de(t)}{dt}\]

Actual Implementation#

There is often reason to change the integral gain \(K_i\) (for example, using different values at different operating points). However, if \(K_i\) changes while \(e_i\) remains unchanged, a large “bump” in \(u(t)\) results. Therefore, the product \(K_i e_i\) is integrated instead.

\[u_i(t) = u_i(t_0) + \int_{t_0}^t K_i(\tau) e(\tau) d\tau\]

The integral error is formulated as a differential equation in the model.

\[\frac{du_i(t)}{dt} = K_i(t) e(t)\]

A naive implementation of derivative action also runs into problems. If the setpoint experiences a step change, the derivative of error is undefined, which can result in numerical problems in simulations and a phenomenon called “derivative kick” in practical implementations. A common way to address this problem is by taking the derivative of \(-y(t)\) rather than \(y_{sp}(t) - y(t)\). The derivative action is the same when the setpoint is fixed but no derivative kick happens when the setpoint is changed. The naive derivative-in-error implementation can be enabled by setting derivative_in_error=True in the config, but it is disabled by default.

Therefore, the final control action for a PID controller with derivative_in_error=False is given by:

\[u(t) = K_p(t) e(t) + u_i(t) - K_d(t) \frac{\text{d}y(t)}{\text{d}t} + u_{ref}(t)\]

Output Limits#

There are two available ways to bound the controller output. The smooth max method provides the closest approximation, but the logistic method provides a method that is smoother and potentially more tractable.

Smooth Max#

Smooth min and smooth max expressions are used to keep the controller output between a minimum and maximum value. The exact form of is given by the following.

\[u_{bound}(t) = \min\left(\max\left(u_{unbound}(t), u_{lb} \right), u_{ub}\right)\]

To avoid a non-smooth formulation, the min and max functions are replaced by smooth versions as follows.

\[\max(a, b) \approx \frac{1}{2}\left[a + b + \left((a-b)^2 + \epsilon^2\right)^{\frac{1}{2}}\right]\]
\[\min(a, b) \approx \frac{1}{2}\left[a + b - \left((a-b)^2 + \epsilon^2\right)^{\frac{1}{2}}\right]\]

Logistic Function#

The logistic bounding function is more smooth and possibly more tractable, but is less faithful to the unbound manipulated variable values and slower to arrive at the bounds for reasonable values of \(k_\text{log}\).

\[u_{bound}(t) = u_{lb} + \frac{(u_{ub} - u_{lb})}{1 + \exp \left( -k_\text{log}\frac{u_{unbound}(t) - (u_{ub} + u_{lb})/2}{(u_{ub} - u_{lb})}\right)}\]

Anti-integral-windup#

When the controller output is clipped to a bound the integral term will continue to grow. This can cause a delay in the controller output moving away from its bound. In order to mitigate this phenomenon, two methods of anti-windup have been implemented: conditional integration and back calculation.

Conditional Integration#

Conditional integration, enabled by setting

antiwindup=ControllerAntiwindupType.CONDITIONAL_INTEGRATION

in the controller config, operates on simple logic:

\[\begin{split}\frac{du_i(t)}{dt} = \begin{cases} K_i(t) e(t) & \text{if} \quad u_{lb}\le u_\text{unbound}\le u_{ub} \\ 0 & \text{otherwise} \end{cases}\end{split}\]

This logic is implemented by first writing this differential equation in terms of Heaviside step functions \(H(x)\).

\[\frac{du_i(t)}{dt} = K_i(t) e(t) \left(H\left(\frac{u_\text{unbound} -u_{lb}}{u_{ub} - u_{lb}}\right) - H\left(\frac{u_\text{unbound} -u_{ub}}{u_{ub} - u_{lb}}\right) \right)\]

and then approximating the (discontinuous) step function using a steep logistic function

\[H(x) \approx \frac{1}{1+\exp(-2k_\text{cond}x)}\]

in which \(k_\text{cond}\) is a parameter governing the steepness of the approximation—the larger it is, the steeper the transitions between 0 and 1 are.

Back Calculation#

Conditional integration is simple to use, but the additional nonlinearity can cause problems for integrators. Back calculation, enabled by setting

antiwindup=ControllerAntiwindupType.BACK_CALCULATION

in the PIDController configuration, does not add any additional nonlinearity, but utilizes the nonlinearity inherent in enforcing controller bounds. The differential equation describing error accumulation is modified:

\[\frac{du_i(t)}{dt} = K_i(t) e(t) + K_b(t)(u_\text{bound}(t) - u_\text{unbound}(t))\]

in which \(K_b(t) \ge 0\) is a new gain term that must be chosen by the user. The more the unbounded control input exceeds the bounded one, the more the integral term unwinds as a result.

To do#

Example with step-by-step initialization#

Dynamic models with a PID controller may be more difficult to initialize. The error and integral error variables may be especially difficult. It is recommended that these models be initialized with a time-stepping method. Examples of this type of initialization should be added.