LIBRARY

PIDEquivalentMPCController

A PID controller in MPC form: without saturation it reproduces DiscreteComponents. DiscretePIDStandard (with its default Forward-Euler integrator and backward-difference filtered derivative) exactly, and its output limits y_min/y_max are handled as input constraints of a QP: whenever any move of the predicted input trajectory would saturate, the control minimizes the accumulated deviation from the PID law over the horizon (given a plant model A, B, C) instead of clamping and back-calculating. Np = 1 reduces exactly to the clamped PID without anti-windup; a longer horizon acts as anti-windup when the saturated spells last long enough to wind the loop up, and makes no difference during brief, deep saturation, where the optimal move is the bound itself (for the control limits the QP itself provides the anti-windup, and the Ni tracking gain — the same parameter and interpretation as in DiscretePIDStandard, but acting on the QP's deviation from the law — extends it to constraints the integrator recursion cannot see, e.g. plant-output constraints).

Construction

The PID with integrator state I and filtered-derivative state D is the static state feedback u = -L x on the augmented state x = [x_p; I; D; r; y] (plant state, controller states, a constant-reference state and the measured output propagated through the plant), with L = [0, -K, -K, -K wp, K]. The MPC minimizes the completed square (u + L x)' R (u + L x), expressed through the weights Q1 = L' R L, Q12 = L' R and Q2 = R of LinearMPCOptimizer — the closed form of inverse-optimal control (Hartley and Maciejowski, Designing MPC controllers by reverse-engineering existing LTI controllers, https://arxiv.org/abs/1109.2816). The optimal cost of the unconstrained problem is identically zero, so the exact LQR cost-to-go (terminal cost) is zero and the equivalence with the PID holds at every prediction horizon Np; the horizon only shapes the constrained behaviour. The integrator and derivative-filter states are computed by the same recursions as DiscretePIDStandard, with the integrator additionally carrying the anti-windup tracking correction (y - v) / (K Ni) in the integrator error — identically zero while no constraint is active, so it never disturbs the equivalence; the plant state is estimated by a steady-state Kalman filter with covariances R1, R2. The current control depends on the measured output, not the plant estimate, so the unconstrained PID equivalence holds regardless of plant-model quality.

Beyond the PID: additional MPC options

The underlying QP accepts more than the PID structure; the following options of LinearMPCController are exposed under the same names (all defaults preserve the exact PID equivalence, and each non-default value trades some of it away, noted below):

• Qy, Qfy: stage and terminal weight on the squared plant-output tracking error (C x_p - r)^2, added to the completed-square PID cost (requires a plant model). Nonzero values bias the controller towards model-based tracking beyond what the PID would have done.

• Q3: weight on the squared control increment.

• Nc: control horizon (the control is held constant after Nc steps); move_block: general move-blocking pattern. The exact equivalence requires the defaults (Nc = Np, no blocking), since blocked moves cannot reproduce the PID trajectory.

• General constraints lb <= Ax x_p + Au u + Ar r <= ub on any linear combination of the plant state, the control and the reference (Ax is m × n_p, Au and Ar are m × 1), with constraint_ks,constraint_soft,constraint_binary,constraint_prioandsoft_weightas inLinearMPCController. Bounds on the plant output are the special caseAx = C: the nameyalready denotes the controller output (the PID convention ofDiscretePIDStandard, where the limitsy_min/y_maxclamp it), so plant-output bounds are written throughAxrather than throughLinearMPCController'symin/ymax.

• binary_controls = [1] with finite y_min/y_max restricts the control to the two limit values (binary / on-off control, relaxable after Nc_binary steps).

• Measured disturbances: Bd (n_p × nd) adds a disturbance input to the plant model and the d input to the component; the model, the constraint predictions and the plant Kalman filter all use it, and disturbance_preview = truetakes the whole known trajectory over the horizon (nd*Np column-major entries) instead of the current value. The unconstrained behaviour remains exactly the PID — the completed-square cost is zero along the PID law whatever the disturbance does — so the channel changes the response only where a constraint is active: the QP then rejects (and, with preview, pre-acts on) the disturbance where a PID can only react through the error. Au_ff-style feed-forward that shifts the unconstrained law itself is not exactly representable in the QP cost forNp > 1 and is intentionally not provided.

Notes

• "Equivalent in the absence of saturation" means the whole predicted input trajectory is unconstrained: near the limits the MPC departs from the PID as soon as any predicted move would saturate, which is what gives it its anti-windup behaviour.

• Without a plant model (A left empty) the component reduces to the clamped PID without anti-windup — supplying at least a rough plant model is recommended whenever the limits can become active.

• All parameters are structural (they are folded into the QP when the model is built) and hence not tunable at runtime. Ts must equal the period of the clock the component is connected to.

• A PID that does not stabilize the supplied plant model is reported with a warning: the unconstrained equivalence still holds, but the constrained behaviour carries no guarantees.

Usage

MPCComponents.PIDEquivalentMPCController(x0=fill(0.0, np_eff), a_d=Td / (Td + Nd * Ts), b_d=Nd * Td / (Td + Nd * Ts), F_p=F, G_p=G, C_p=Cp, Gd_p=Gd, K_kf=_pidmpc_kf_K(F, G, Cp, np, R1, R2), d_dummy=fill(0.0, nd_in_eff))

Parameters:

Name	Description	Units	Default value
A	Plant state matrix (n_p × n_p); empty for no plant model (clamped PID without anti-windup)	–	fill(0.0, 0, 0)
B	Plant input matrix (n_p × 1)	–	fill(0.0, 0, 0)
C	Plant output (measurement) matrix (1 × n_p); empty for the identity when n_p = 1	–	fill(0.0, 0, 0)
Bd	Plant input matrix of the measured disturbances (n_p × nd); empty for none	–	fill(0.0, 0, 0)
continuous	Interpret A and B as continuous-time and ZOH-discretize them with Ts	–	false
Ts	Sample time; must equal the period of the clock the component is connected to	–	1.0
Np	Prediction horizon (only affects the behaviour when the output limits become active)	–	10
Nc	Control horizon; the control is held constant after Nc steps. The exact PID equivalence requires the default Nc = Np	–	Np
disturbance_preview	Take a measured-disturbance trajectory over the horizon (nd*Np entries) instead of the current value	–	false
with_I	Whether or not to include the integral part	–	true
with_D	Whether or not to include the derivative part	–	true
K	Proportional gain	–	1.0
Ti	Integral time constant (only active if with_I = true)	–	1.0
Td	Derivative time (only active if with_D = true)	–	1.0
Nd	Maximum derivative gain (only active if with_D = true). Typically set to 10-100.	–	10.0
wp	Set-point weighting in the proportional part. Set to 0 to prevent step changes in the output due to step changes in the reference.	–	1.0
wd	Set-point weighting in the derivative part. Set to 0 to prevent very large impulsive changes in the output due to step changes in the reference.	–	1.0
Ni	Anti-windup tracking gain, as in DiscretePIDStandard: the integrator error becomes eI = e + (y - v) / (K Ni), where y is the applied control and v the unconstrained PID law, so the correction is zero whenever no constraint is active (the exact PID equivalence is untouched) and drains the integrator whenever the QP deviates from the law. Inf (the default, unlike DiscretePIDStandard's sqrt(Td/Ti)) disables it: enabling it is recommended whenever output/state constraints are used — without it the integrator winds up against them — but keep it disabled when cost-modifying options (Qy, Qfy, Q3, Nc < Np, move_block) are active, since those deviate from the PID law by design and the correction would fight them. Deliberately kept out of the QP's prediction model, which would let the optimizer steer the integrator and relax its own target law	–	Inf
y_max	Maximum controller output	–	Inf
y_min	Minimum controller output	–	ifelse(y_ma..._max, -Inf)
R_cost	Weight of the completed-square cost; any positive value gives the same unconstrained controller, the value scales the QP	–	1.0
Qy	Stage weight on the squared plant-output tracking error (C x_p - r)^2, added to the PID cost; 0 preserves the exact PID equivalence (requires a plant model)	–	0.0
Qfy	Terminal weight on the squared plant-output tracking error; 0 preserves the exact PID equivalence (requires a plant model)	–	0.0
Q3	Weight on the squared control increment; 0 preserves the exact PID equivalence	–	0.0
Ax	Plant-state coefficient of the general constraints lb <= Ax x_p + Au u + Ar r <= ub (m × n_p); empty for zero. Plant-output bounds are the special case Ax = C	–	fill(0.0, 0, 0)
Au	Control coefficient of the general constraints (m × 1); empty for zero	–	fill(0.0, 0, 0)
Ar	Reference coefficient of the general constraints (m × 1); empty for zero	–	fill(0.0, 0, 0)
lb	Lower bounds of the general constraints (m); empty for none	–	fill(0.0, 0)
ub	Upper bounds of the general constraints (m); empty for none	–	fill(0.0, 0)
constraint_ks	Prediction steps on which the general constraints are enforced; empty for 2:Np	–	[]
constraint_soft	Soften the general constraints	–	false
constraint_binary	Enforce the general constraints as binary (one of the bounds holds with equality)	–	false
constraint_prio	Priority of the general constraints	–	0
soft_weight	Penalty weight on soft-constraint violations	–	1e6
primal_tol	DAQP primal feasibility tolerance; -1 keeps the solver default. Applied to the Julia solver and embedded in the generated C	–	-1.0
dual_tol	DAQP dual feasibility (optimality) tolerance; -1 keeps the solver default	–	-1.0
iter_limit	DAQP iteration limit; -1 keeps the solver default	–	-1
binary_controls	Indices of the controls restricted to the values {y_min, y_max}; empty for none. The only valid index is 1	–	[]
Nc_binary	Number of initial horizon steps over which the binary controls are restricted to {y_min, y_max}, the rest are relaxed to continuous; -1 for the whole control horizon	–	-1
move_block	Move-blocking pattern, the number of steps each control move is held; empty for none. Non-default blocking breaks the exact PID equivalence	–	[]
R1	Process-noise covariance of the plant Kalman filter (n_p × n_p, or n_p diagonal entries); empty for the identity	–	fill(0.0, 0, 0)
R2	Measurement-noise covariance of the plant Kalman filter (1 × 1); empty for the identity	–	fill(0.0, 0, 0)
backend	QP solver backend; Auto uses the compiled C code when a C compiler is available	–	Backend.Auto()
precision	Float type of the generated QP code (C backend only)	–	MPCPrecision.Double()
output_trajectories	Also record the predicted state and control trajectories and the QP statistics of every tick (variables optimizer.x_pred, optimizer.u_pred, optimizer.stats), for debugging and mpc_gui; Julia backend only. The readbacks are three equations whose number does not grow with the prediction horizon, though the trajectories themselves are stored in the solution; output_statistics records the statistics alone	–	false
output_statistics	Record the QP statistics of every tick (variable optimizer.stats) without the predicted trajectories; Julia backend only. One equation. Implied by output_trajectories	–	false
np	Plant state dimension (derived)	–	size(A, 1)
nd	Number of measured disturbances (derived)	–	size(Bd, 2)
has_dist	Whether measured disturbances are present (derived)	–	nd > 0
nd_in	Width of the disturbance input: nd, or nd*Np with disturbance preview (derived)	–	ifelse(dist...d * Np, nd)
nd_in_eff	Width of the internal disturbance vector, max(nd_in, 1) (derived)	–	max(nd_in, 1)
has_plant	Whether a plant model was given (derived)	–	np > 0
np_eff	Width of the plant-estimator vectors, max(np, 1) (derived)	–	max(np, 1)
nxt	State dimension of the augmented model [x_p; I; D; r; y] (derived)	–	np + ifelse..., 1, 0) + 2
i_I	Index of the integrator state I in the augmented state (derived)	–	np + 1
i_D	Index of the derivative-filter state D in the augmented state (derived)	–	np + ifelse..., 1, 0) + 1
i_r	Index of the constant-reference state r in the augmented state (derived)	–	nxt - 1
i_y	Index of the measured-output state y in the augmented state (derived)	–	nxt
F	Discrete-time plant state matrix, padded to np_eff × np_eff (derived)	–	_pidmpc_F(A...inuous, Ts)
G	Discrete-time plant input matrix, padded to np_eff × 1 (derived)	–	_pidmpc_G(A...inuous, Ts)
Cp	Plant output matrix, padded to 1 × np_eff (derived)	–	_pidmpc_C(A...inuous, Ts)
Gd	Discrete-time plant disturbance input matrix, padded to np_eff × nd_in_eff; only the first nd columns (the current disturbance) are nonzero (derived)	–	_pidmpc_Gd(... nd_in_eff)
Atilde	Augmented model state matrix (derived)	–	_pidmpc_A(F..._I, with_D)
Btilde	Augmented model input matrix (derived)	–	_pidmpc_B(G..._I, with_D)
Btilde_d	Augmented model disturbance input matrix (derived)	–	_pidmpc_Bd(..._I, with_D)
Ltilde	The PID law as state feedback u = -Ltilde x, Ltilde = [0, -K, -K, -K wp, K] (derived)	–	_pidmpc_L(A..._I, with_D)
Ttilde	Row extracting the plant-output tracking error C x_p - r from the augmented state (derived)	–	_pidmpc_tra...D, Qy, Qfy)
Q1t	Stage state weight of the QP: the completed square R_cost (u + Ltilde x)^2 contributes R_cost Ltilde' Ltilde, the optional plant-output tracking cost Qy (C x_p - r)^2 adds Qy Ttilde' Ttilde (derived)	–	R_cost * tr...e) * Ttilde
Q12t	State-control cross weight of the completed square, R_cost Ltilde' (derived)	–	R_cost * tr...ose(Ltilde)
Q2t	Control weight of the completed square, R_cost (derived)	–	fill(R_cost, 1, 1)
Qfxt	Terminal state weight, the plant-output tracking cost Qfy (C x_p - r)^2 at the end of the horizon (derived)	–	Qfy * trans...e) * Ttilde
mpc	LinearMPC.jl controller with the completed-square cost (derived)	–	build_linea...bels=["y"])
x0	Initial plant state estimate	–	fill(0.0, np_eff)

Connectors

• u_s - This connector represents a real signal as an input to a component (RealInput)

• u_m - This connector represents a real signal as an input to a component (RealInput)

• d - This connector represents a real signal as an input to a component (RealInput)

• y - This connector represents a real signal as an output from a component (RealOutput)

• exitflag - This connector represents a real signal as an output from a component (RealOutput)

Variables

Name	Description	Units
e	Error signal	–
wde	Setpoint-weighted error for the derivative	–
I_state	State of the integrator	–
D_state	State of the filtered derivative	–
xcorr	Corrected plant state estimate	–
xpred	Predicted plant state estimate	–
ym	Measurement as a vector, for the estimator equations	–
uv	Applied control as a vector, for the estimator equations	–
d_int	Measured disturbances as a vector, for the estimator equations	–
v_pid	Unconstrained PID law, for the anti-windup tracking correction	–
eI	Integrator error including the anti-windup tracking signal	–

Behavior

Behavior of this component cannot be rendered because it includes path variables.

Source

"""
A PID controller in MPC form: without saturation it reproduces `DiscreteComponents.
DiscretePIDStandard` (with its default Forward-Euler integrator and backward-difference filtered
derivative) exactly, and its output limits `y_min`/`y_max` are handled as input constraints of a
QP: whenever any move of the *predicted* input trajectory would saturate, the control minimizes
the accumulated deviation from the PID law over the horizon (given a plant model `A`, `B`, `C`)
instead of clamping and back-calculating. `Np = 1` reduces exactly to the clamped PID without
anti-windup; a longer horizon acts as anti-windup when the saturated spells last long enough to
wind the loop up, and makes no difference during brief, deep saturation, where the optimal move
is the bound itself (for the control
limits the QP itself provides the anti-windup, and the `Ni` tracking gain — the same parameter
and interpretation as in `DiscretePIDStandard`, but acting on the QP's deviation from the law —
extends it to constraints the integrator recursion cannot see, e.g. plant-output constraints).

# Construction
The PID with integrator state `I` and filtered-derivative state `D` is the static state feedback
`u = -L x` on the augmented state `x = [x_p; I; D; r; y]` (plant state, controller states, a
constant-reference state and the measured output propagated through the plant), with
`L = [0, -K, -K, -K wp, K]`. The MPC minimizes the completed square `(u + L x)' R (u + L x)`,
expressed through the weights `Q1 = L' R L`, `Q12 = L' R` and `Q2 = R` of `LinearMPCOptimizer` —
the closed form of inverse-optimal control (Hartley and Maciejowski, *Designing MPC controllers by
reverse-engineering existing LTI controllers*, https://arxiv.org/abs/1109.2816). The optimal cost
of the unconstrained problem is identically zero, so the exact LQR cost-to-go (terminal cost) is
zero and the equivalence with the PID holds at **every** prediction horizon `Np`; the horizon only
shapes the constrained behaviour. The integrator and derivative-filter states are computed
by the same recursions as `DiscretePIDStandard`, with the integrator additionally carrying the
anti-windup tracking correction `(y - v) / (K Ni)` in the integrator error — identically zero while no constraint is
active, so it never disturbs the equivalence; the plant state is estimated by a steady-state
Kalman filter with covariances `R1`, `R2`. The current control depends on the *measured* output,
not the plant estimate, so the unconstrained PID equivalence holds regardless of plant-model
quality.

# Beyond the PID: additional MPC options
The underlying QP accepts more than the PID structure; the following options of
`LinearMPCController` are exposed under the same names (all defaults preserve the exact PID
equivalence, and each non-default value trades some of it away, noted below):
- `Qy`, `Qfy`: stage and terminal weight on the squared plant-output tracking error
  `(C x_p - r)^2`, added to the completed-square PID cost (requires a plant model). Nonzero values
  bias the controller towards model-based tracking beyond what the PID would have done.
- `Q3`: weight on the squared control increment.
- `Nc`: control horizon (the control is held constant after `Nc` steps); `move_block`: general
  move-blocking pattern. The exact equivalence requires the defaults (`Nc = Np`, no blocking),
  since blocked moves cannot reproduce the PID trajectory.
- General constraints `lb <= Ax x_p + Au u + Ar r <= ub` on any linear combination of the plant
  state, the control and the reference (`Ax` is m × n_p, `Au` and `Ar` are m × 1), with
  `constraint_ks`, `constraint_soft`, `constraint_binary`, `constraint_prio` and `soft_weight` as
  in `LinearMPCController`. Bounds on the *plant output* are the special case `Ax = C`: the name
  `y` already denotes the controller output (the PID convention of `DiscretePIDStandard`, where
  the limits `y_min`/`y_max` clamp it), so plant-output bounds are written through `Ax` rather
  than through `LinearMPCController`'s `ymin`/`ymax`.
- `binary_controls = [1]` with finite `y_min`/`y_max` restricts the control to the two limit
  values (binary / on-off control, relaxable after `Nc_binary` steps).
- Measured disturbances: `Bd` (n_p × nd) adds a disturbance input to the plant model and the `d`
  input to the component; the model, the constraint predictions and the plant Kalman filter all
  use it, and `disturbance_preview = true` takes the whole known trajectory over the horizon
  (nd*Np column-major entries) instead of the current value. The unconstrained behaviour remains
  *exactly* the PID — the completed-square cost is zero along the PID law whatever the
  disturbance does — so the channel changes the response only where a constraint is active: the
  QP then rejects (and, with preview, pre-acts on) the disturbance where a PID can only react
  through the error. A `u_ff`-style feed-forward that shifts the unconstrained law itself is not
  exactly representable in the QP cost for `Np > 1` and is intentionally not provided.

# Notes
- "Equivalent in the absence of saturation" means the whole *predicted* input trajectory is
  unconstrained: near the limits the MPC departs from the PID as soon as any predicted move would
  saturate, which is what gives it its anti-windup behaviour.
- Without a plant model (`A` left empty) the component reduces to the clamped PID without
  anti-windup — supplying at least a rough plant model is recommended whenever the limits can
  become active.
- All parameters are structural (they are folded into the QP when the model is built) and hence
  not tunable at runtime. `Ts` must equal the period of the clock the component is connected to.
- A PID that does not stabilize the supplied plant model is reported with a warning: the
  unconstrained equivalence still holds, but the constrained behaviour carries no guarantees.
"""
component PIDEquivalentMPCController@[input clk extends Discrete]
  "The reference signal to the controller"
  u_s = RealInput@[clk]() {
    "Dyad": {"placement": {"diagram": {"x1": -100, "y1": 200, "x2": 0, "y2": 300}}}
  }
  "The measurement feedback"
  u_m = RealInput@[clk]() {
    "Dyad": {"placement": {"diagram": {"x1": -100, "y1": 700, "x2": 0, "y2": 800}}}
  }
  "Measured disturbances; nd entries, or nd*Np column-major entries with disturbance preview (only when Bd is given)"
  d = [RealInput@[clk]() for i in 1:nd_in] if has_dist {
    "Dyad": {
      "placement": {"diagram": {"x1": 450, "y1": -100, "x2": 550, "y2": 0, "rot": 90}}
    }
  }
  "The control signal output"
  y = RealOutput@[clk]() {
    "Dyad": {"placement": {"diagram": {"x1": 1000, "y1": 700, "x2": 1100, "y2": 800}}}
  }
  "Solver status of the last step (DAQP exit flag: 1 optimal, 2 with active soft constraints, negative on failure); lmpc_status_name names it"
  exitflag = RealOutput@[clk]() {
    "Dyad": {"placement": {"diagram": {"x1": 1000, "y1": 200, "x2": 1100, "y2": 300}}}
  }
  "Constrained optimizer solving the completed-square QP"
  optimizer = MPCComponents.LinearMPCOptimizer@[clk = clk](mpc = mpc, backend = backend, precision = precision, output_trajectories = output_trajectories, output_statistics = output_statistics) {
    "Dyad": {"placement": {"diagram": {"x1": 400, "y1": 400, "x2": 600, "y2": 600}}}
  }
  "Plant state matrix (n_p × n_p); empty for no plant model (clamped PID without anti-windup)"
  structural parameter A::Real[:, :] = fill(0.0, 0, 0)
  "Plant input matrix (n_p × 1)"
  structural parameter B::Real[:, :] = fill(0.0, 0, 0)
  "Plant output (measurement) matrix (1 × n_p); empty for the identity when n_p = 1"
  structural parameter C::Real[:, :] = fill(0.0, 0, 0)
  "Plant input matrix of the measured disturbances (n_p × nd); empty for none"
  structural parameter Bd::Real[:, :] = fill(0.0, 0, 0)
  "Interpret A and B as continuous-time and ZOH-discretize them with Ts"
  structural parameter continuous::Boolean = false
  "Sample time; must equal the period of the clock the component is connected to"
  structural parameter Ts::Real = 1.0
  "Prediction horizon (only affects the behaviour when the output limits become active)"
  structural parameter Np::Integer = 10
  "Control horizon; the control is held constant after Nc steps. The exact PID equivalence requires the default Nc = Np"
  structural parameter Nc::Integer = Np
  "Take a measured-disturbance trajectory over the horizon (nd*Np entries) instead of the current value"
  structural parameter disturbance_preview::Boolean = false
  "Whether or not to include the integral part"
  structural parameter with_I::Boolean = true
  "Whether or not to include the derivative part"
  structural parameter with_D::Boolean = true
  "Proportional gain"
  structural parameter K::Real = 1.0
  "Integral time constant (only active if `with_I = true`)"
  structural parameter Ti::Real = 1.0
  "Derivative time (only active if `with_D = true`)"
  structural parameter Td::Real = 1.0
  "Maximum derivative gain (only active if `with_D = true`). Typically set to 10-100."
  structural parameter Nd::Real = 10.0
  "Set-point weighting in the proportional part. Set to `0` to prevent step changes in the output due to step changes in the reference."
  structural parameter wp::Real = 1.0
  "Set-point weighting in the derivative part. Set to `0` to prevent very large impulsive changes in the output due to step changes in the reference."
  structural parameter wd::Real = 1.0
  "Anti-windup tracking gain, as in DiscretePIDStandard: the integrator error becomes eI = e + (y - v) / (K Ni), where y is the applied control and v the unconstrained PID law, so the correction is zero whenever no constraint is active (the exact PID equivalence is untouched) and drains the integrator whenever the QP deviates from the law. Inf (the default, unlike DiscretePIDStandard's sqrt(Td/Ti)) disables it: enabling it is recommended whenever output/state constraints are used — without it the integrator winds up against them — but keep it disabled when cost-modifying options (Qy, Qfy, Q3, Nc < Np, move_block) are active, since those deviate from the PID law by design and the correction would fight them. Deliberately kept out of the QP's prediction model, which would let the optimizer steer the integrator and relax its own target law"
  structural parameter Ni::Real = Inf
  "Maximum controller output"
  structural parameter y_max::Real = Inf
  "Minimum controller output"
  structural parameter y_min::Real = ifelse(y_max > 0, -y_max, -Inf)
  "Weight of the completed-square cost; any positive value gives the same unconstrained controller, the value scales the QP"
  structural parameter R_cost::Real = 1.0
  "Stage weight on the squared plant-output tracking error (C x_p - r)^2, added to the PID cost; 0 preserves the exact PID equivalence (requires a plant model)"
  structural parameter Qy::Real = 0.0
  "Terminal weight on the squared plant-output tracking error; 0 preserves the exact PID equivalence (requires a plant model)"
  structural parameter Qfy::Real = 0.0
  "Weight on the squared control increment; 0 preserves the exact PID equivalence"
  structural parameter Q3::Real = 0.0
  "Plant-state coefficient of the general constraints lb <= Ax x_p + Au u + Ar r <= ub (m × n_p); empty for zero. Plant-output bounds are the special case Ax = C"
  structural parameter Ax::Real[:, :] = fill(0.0, 0, 0)
  "Control coefficient of the general constraints (m × 1); empty for zero"
  structural parameter Au::Real[:, :] = fill(0.0, 0, 0)
  "Reference coefficient of the general constraints (m × 1); empty for zero"
  structural parameter Ar::Real[:, :] = fill(0.0, 0, 0)
  "Lower bounds of the general constraints (m); empty for none"
  structural parameter lb::Real[:] = fill(0.0, 0)
  "Upper bounds of the general constraints (m); empty for none"
  structural parameter ub::Real[:] = fill(0.0, 0)
  "Prediction steps on which the general constraints are enforced; empty for 2:Np"
  structural parameter constraint_ks::Integer[:] = []
  "Soften the general constraints"
  structural parameter constraint_soft::Boolean = false
  "Enforce the general constraints as binary (one of the bounds holds with equality)"
  structural parameter constraint_binary::Boolean = false
  "Priority of the general constraints"
  structural parameter constraint_prio::Integer = 0
  "Penalty weight on soft-constraint violations"
  structural parameter soft_weight::Real = 1e6
  "DAQP primal feasibility tolerance; -1 keeps the solver default. Applied to the Julia solver and embedded in the generated C"
  structural parameter primal_tol::Real = -1.0
  "DAQP dual feasibility (optimality) tolerance; -1 keeps the solver default"
  structural parameter dual_tol::Real = -1.0
  "DAQP iteration limit; -1 keeps the solver default"
  structural parameter iter_limit::Integer = -1
  "Indices of the controls restricted to the values {y_min, y_max}; empty for none. The only valid index is 1"
  structural parameter binary_controls::Integer[:] = []
  "Number of initial horizon steps over which the binary controls are restricted to {y_min, y_max}, the rest are relaxed to continuous; -1 for the whole control horizon"
  structural parameter Nc_binary::Integer = -1
  "Move-blocking pattern, the number of steps each control move is held; empty for none. Non-default blocking breaks the exact PID equivalence"
  structural parameter move_block::Integer[:] = []
  "Process-noise covariance of the plant Kalman filter (n_p × n_p, or n_p diagonal entries); empty for the identity"
  structural parameter R1::Real[:, :] = fill(0.0, 0, 0)
  "Measurement-noise covariance of the plant Kalman filter (1 × 1); empty for the identity"
  structural parameter R2::Real[:, :] = fill(0.0, 0, 0)
  "QP solver backend; Auto uses the compiled C code when a C compiler is available"
  structural parameter backend::Backend = Backend.Auto()
  "Float type of the generated QP code (C backend only)"
  structural parameter precision::MPCPrecision = MPCPrecision.Double()
  "Also record the predicted state and control trajectories and the QP statistics of every tick (variables optimizer.x_pred, optimizer.u_pred, optimizer.stats), for debugging and mpc_gui; Julia backend only. The readbacks are three equations whose number does not grow with the prediction horizon, though the trajectories themselves are stored in the solution; output_statistics records the statistics alone"
  structural parameter output_trajectories::Boolean = false
  "Record the QP statistics of every tick (variable optimizer.stats) without the predicted trajectories; Julia backend only. One equation. Implied by output_trajectories"
  structural parameter output_statistics::Boolean = false
  "Plant state dimension (derived)"
  final structural parameter np::Integer = size(A, 1)
  "Number of measured disturbances (derived)"
  final structural parameter nd::Integer = size(Bd, 2)
  "Whether measured disturbances are present (derived)"
  final structural parameter has_dist::Boolean = nd > 0
  "Width of the disturbance input: nd, or nd*Np with disturbance preview (derived)"
  final structural parameter nd_in::Integer = ifelse(disturbance_preview, nd * Np, nd)
  "Width of the internal disturbance vector, max(nd_in, 1) (derived)"
  final structural parameter nd_in_eff::Integer = max(nd_in, 1)
  "Whether a plant model was given (derived)"
  final structural parameter has_plant::Boolean = np > 0
  "Width of the plant-estimator vectors, max(np, 1) (derived)"
  final structural parameter np_eff::Integer = max(np, 1)
  "State dimension of the augmented model [x_p; I; D; r; y] (derived)"
  final structural parameter nxt::Integer = np + ifelse(with_I, 1, 0) + ifelse(with_D, 1, 0) + 2
  "Index of the integrator state I in the augmented state (derived)"
  final structural parameter i_I::Integer = np + 1
  "Index of the derivative-filter state D in the augmented state (derived)"
  final structural parameter i_D::Integer = np + ifelse(with_I, 1, 0) + 1
  "Index of the constant-reference state r in the augmented state (derived)"
  final structural parameter i_r::Integer = nxt - 1
  "Index of the measured-output state y in the augmented state (derived)"
  final structural parameter i_y::Integer = nxt
  "Discrete-time plant state matrix, padded to np_eff × np_eff (derived)"
  final structural parameter F::Real[np_eff, np_eff] = _pidmpc_F(A, B, C, continuous, Ts)
  "Discrete-time plant input matrix, padded to np_eff × 1 (derived)"
  final structural parameter G::Real[np_eff, 1] = _pidmpc_G(A, B, C, continuous, Ts)
  "Plant output matrix, padded to 1 × np_eff (derived)"
  final structural parameter Cp::Real[1, np_eff] = _pidmpc_C(A, B, C, continuous, Ts)
  "Discrete-time plant disturbance input matrix, padded to np_eff × nd_in_eff; only the first nd columns (the current disturbance) are nonzero (derived)"
  final structural parameter Gd::Real[np_eff, nd_in_eff] = _pidmpc_Gd(A, B, C, Bd, continuous, Ts, nd_in_eff)
  "Augmented model state matrix (derived)"
  final structural parameter Atilde::Real[nxt, nxt] = _pidmpc_A(F, G, Cp, np, Ts, Ti, Td, Nd, with_I, with_D)
  "Augmented model input matrix (derived)"
  final structural parameter Btilde::Real[nxt, 1] = _pidmpc_B(G, Cp, np, Ts, Td, Nd, with_I, with_D)
  "Augmented model disturbance input matrix (derived)"
  final structural parameter Btilde_d::Real[nxt, nd] = _pidmpc_Bd(Gd, Cp, np, nd, Ts, Td, Nd, with_I, with_D)
  "The PID law as state feedback u = -Ltilde x, Ltilde = [0, -K, -K, -K wp, K] (derived)"
  final structural parameter Ltilde::Real[1, nxt] = _pidmpc_L(Atilde, Btilde, np, K, wp, with_I, with_D)
  "Row extracting the plant-output tracking error C x_p - r from the augmented state (derived)"
  final structural parameter Ttilde::Real[1, nxt] = _pidmpc_tracking_row(Cp, np, with_I, with_D, Qy, Qfy)
  "Stage state weight of the QP: the completed square R_cost (u + Ltilde x)^2 contributes R_cost Ltilde' Ltilde, the optional plant-output tracking cost Qy (C x_p - r)^2 adds Qy Ttilde' Ttilde (derived)"
  final structural parameter Q1t::Real[nxt, nxt] = R_cost * transpose(Ltilde) * Ltilde + Qy * transpose(Ttilde) * Ttilde
  "State-control cross weight of the completed square, R_cost Ltilde' (derived)"
  final structural parameter Q12t::Real[nxt, 1] = R_cost * transpose(Ltilde)
  "Control weight of the completed square, R_cost (derived)"
  final structural parameter Q2t::Real[1, 1] = fill(R_cost, 1, 1)
  "Terminal state weight, the plant-output tracking cost Qfy (C x_p - r)^2 at the end of the horizon (derived)"
  final structural parameter Qfxt::Real[nxt, nxt] = Qfy * transpose(Ttilde) * Ttilde
  "LinearMPC.jl controller with the completed-square cost (derived)"
  final structural parameter mpc::MPCRef = build_linear_mpc(A = Atilde, B = Btilde, Bd = Btilde_d, disturbance_preview = disturbance_preview, Ts = Ts, Np = Np, Nc = Nc, Q1 = Q1t, Q2 = Q2t, Q3 = _pidmpc_Q3(Q3), Q12 = Q12t, Qf = _pidmpc_Qf(nxt), Qfx = Qfxt, reference_tracking = false, umin = _pidmpc_umin(y_min, y_max), umax = _pidmpc_umax(y_min, y_max), Ax = _pidmpc_Ax(Ax, Ar, np, with_I, with_D), Au = Au, lb = lb, ub = ub, constraint_ks = constraint_ks, constraint_soft = constraint_soft, constraint_binary = constraint_binary, constraint_prio = constraint_prio, soft_weight = soft_weight, primal_tol = primal_tol, dual_tol = dual_tol, iter_limit = iter_limit, binary_controls = binary_controls, Nc_binary = Nc_binary, move_block = move_block, state_labels = _pidmpc_state_labels(np, with_I, with_D), control_labels = ["y"])
  "Initial plant state estimate"
  parameter x0::Real[np_eff] = fill(0.0, np_eff)
  "Derivative-filter pole, as in the backward-difference filter of DiscretePIDStandard (derived)"
  final parameter a_d::Real = Td / (Td + Nd * Ts)
  "Derivative-filter gain (derived)"
  final parameter b_d::Real = Nd * Td / (Td + Nd * Ts)
  "Plant estimator state matrix, F as a runtime parameter (derived)"
  final parameter F_p::Real[np_eff, np_eff] = F
  "Plant estimator input matrix, G as a runtime parameter (derived)"
  final parameter G_p::Real[np_eff, 1] = G
  "Plant estimator measurement matrix, Cp as a runtime parameter (derived)"
  final parameter C_p::Real[1, np_eff] = Cp
  "Plant estimator disturbance input matrix, Gd as a runtime parameter (derived)"
  final parameter Gd_p::Real[np_eff, nd_in_eff] = Gd
  "Steady-state Kalman gain of the plant estimator (derived)"
  final parameter K_kf::Real[np_eff, 1] = _pidmpc_kf_K(F, G, Cp, np, R1, R2)
  "Zero placeholder for the disturbance vector when d is absent"
  final parameter d_dummy::Real[nd_in_eff] = fill(0.0, nd_in_eff)
  "Error signal"
  variable e::Real
  "Setpoint-weighted error for the derivative"
  variable wde::Real
  "State of the integrator"
  variable I_state::Real
  "State of the filtered derivative"
  variable D_state::Real
  "Corrected plant state estimate"
  variable xcorr::Real[np_eff]
  "Predicted plant state estimate"
  variable xpred::Real[np_eff]
  "Measurement as a vector, for the estimator equations"
  variable ym::Real[1]
  "Applied control as a vector, for the estimator equations"
  variable uv::Real[1]
  "Measured disturbances as a vector, for the estimator equations"
  variable d_int::Real[nd_in_eff]
  "Unconstrained PID law, for the anti-windup tracking correction"
  variable v_pid::Real
  "Integrator error including the anti-windup tracking signal"
  variable eI::Real
relations
  e = u_s - u_m
  wde = wd * u_s - u_m
  v_pid = K * (wp * u_s - u_m + I_state + D_state)
  eI = e + (y - v_pid) / (K * Ni)
  ym[1] = u_m
  uv[1] = y
  # The augmented state [x_p; I; D; r; y] handed to the optimizer
  if with_I
    initial I_state@(clk-1) = 0
    initial eI@(clk-1) = 0
    I_state@clk = I_state@(clk-1) + Ts / Ti * eI@(clk-1)
    optimizer.xhat[i_I] = I_state
  else
    I_state = 0
  end
  if with_D
    initial D_state@(clk-1) = 0
    initial wde@(clk-1) = 0
    D_state@clk = a_d * D_state@(clk-1) + b_d * (wde@clk - wde@(clk-1))
    optimizer.xhat[i_D] = D_state
  else
    D_state = 0
  end
  optimizer.xhat[i_r] = u_s
  optimizer.xhat[i_y] = u_m
  if has_dist
    d_int = d
  else
    d_int = d_dummy
  end
  if has_plant
    initial xpred@(clk-1) = x0
    xcorr@clk = xpred@(clk-1) + K_kf * (ym@clk - C_p * xpred@(clk-1))
    xpred@clk = F_p * xcorr@clk + G_p * uv@clk + Gd_p * d_int@clk
    for i in 1:np
      optimizer.xhat[i] = xcorr[i]
    end
  else
    xcorr = fill(0.0, np_eff)
    xpred = fill(0.0, np_eff)
  end
  if has_dist
    for i in 1:nd_in
      connect(d[i], optimizer.d[i])
    end
  end
  connect(optimizer.u[1], y) {}
  connect(optimizer.exitflag, exitflag) {}
metadata {
  "Dyad": {
    "icons": {"default": "dyad://MPCComponents/PIDEquivalentMPCController.svg"},
    "doc": {"behavior": false},
    "panel": {
      "control": "tabs",
      "contents": [
        "K",
        "Ti",
        "Td",
        "Nd",
        "wp",
        "wd",
        "Ni",
        "y_max",
        "y_min",
        "with_I",
        "with_D",
        "Ts"
      ],
      "children": {
        "Plant Model": {"contents": ["A", "B", "C", "Bd", "continuous"]},
        "State Estimator": {"contents": ["R1", "R2", "x0"]},
        "MPC": {
          "children": {
            "Horizon": {"contents": ["Np", "Nc", "move_block", "disturbance_preview"]},
            "Cost": {"contents": ["R_cost", "Qy", "Qfy", "Q3"]},
            "Constraints": {
              "contents": [
                "Ax",
                "Au",
                "Ar",
                "lb",
                "ub",
                "constraint_ks",
                "constraint_soft",
                "constraint_binary",
                "constraint_prio",
                "soft_weight",
                "binary_controls",
                "Nc_binary"
              ]
            }
          }
        },
        "Solver": {
          "contents": [
            "backend",
            "precision",
            "output_trajectories",
            "output_statistics",
            "primal_tol",
            "dual_tol",
            "iter_limit"
          ]
        },
        "Internal": {
          "contents": [
            "np",
            "nd",
            "has_dist",
            "nd_in",
            "nd_in_eff",
            "has_plant",
            "np_eff",
            "nxt",
            "i_I",
            "i_D",
            "i_r",
            "i_y",
            "F",
            "G",
            "Cp",
            "Gd",
            "Btilde_d",
            "d_dummy",
            "d_int",
            "Atilde",
            "Btilde",
            "Ltilde",
            "Ttilde",
            "Q1t",
            "Q12t",
            "Q2t",
            "Qfxt",
            "mpc",
            "a_d",
            "b_d",
            "v_pid",
            "eI",
            "F_p",
            "G_p",
            "C_p",
            "Gd_p",
            "K_kf",
            "e",
            "wde",
            "I_state",
            "D_state",
            "xcorr",
            "xpred",
            "ym",
            "uv"
          ]
        }
      }
    }
  }
}
end

Flattened Source

"""
A PID controller in MPC form: without saturation it reproduces `DiscreteComponents.
DiscretePIDStandard` (with its default Forward-Euler integrator and backward-difference filtered
derivative) exactly, and its output limits `y_min`/`y_max` are handled as input constraints of a
QP: whenever any move of the *predicted* input trajectory would saturate, the control minimizes
the accumulated deviation from the PID law over the horizon (given a plant model `A`, `B`, `C`)
instead of clamping and back-calculating. `Np = 1` reduces exactly to the clamped PID without
anti-windup; a longer horizon acts as anti-windup when the saturated spells last long enough to
wind the loop up, and makes no difference during brief, deep saturation, where the optimal move
is the bound itself (for the control
limits the QP itself provides the anti-windup, and the `Ni` tracking gain — the same parameter
and interpretation as in `DiscretePIDStandard`, but acting on the QP's deviation from the law —
extends it to constraints the integrator recursion cannot see, e.g. plant-output constraints).

# Construction
The PID with integrator state `I` and filtered-derivative state `D` is the static state feedback
`u = -L x` on the augmented state `x = [x_p; I; D; r; y]` (plant state, controller states, a
constant-reference state and the measured output propagated through the plant), with
`L = [0, -K, -K, -K wp, K]`. The MPC minimizes the completed square `(u + L x)' R (u + L x)`,
expressed through the weights `Q1 = L' R L`, `Q12 = L' R` and `Q2 = R` of `LinearMPCOptimizer` —
the closed form of inverse-optimal control (Hartley and Maciejowski, *Designing MPC controllers by
reverse-engineering existing LTI controllers*, https://arxiv.org/abs/1109.2816). The optimal cost
of the unconstrained problem is identically zero, so the exact LQR cost-to-go (terminal cost) is
zero and the equivalence with the PID holds at **every** prediction horizon `Np`; the horizon only
shapes the constrained behaviour. The integrator and derivative-filter states are computed
by the same recursions as `DiscretePIDStandard`, with the integrator additionally carrying the
anti-windup tracking correction `(y - v) / (K Ni)` in the integrator error — identically zero while no constraint is
active, so it never disturbs the equivalence; the plant state is estimated by a steady-state
Kalman filter with covariances `R1`, `R2`. The current control depends on the *measured* output,
not the plant estimate, so the unconstrained PID equivalence holds regardless of plant-model
quality.

# Beyond the PID: additional MPC options
The underlying QP accepts more than the PID structure; the following options of
`LinearMPCController` are exposed under the same names (all defaults preserve the exact PID
equivalence, and each non-default value trades some of it away, noted below):
- `Qy`, `Qfy`: stage and terminal weight on the squared plant-output tracking error
  `(C x_p - r)^2`, added to the completed-square PID cost (requires a plant model). Nonzero values
  bias the controller towards model-based tracking beyond what the PID would have done.
- `Q3`: weight on the squared control increment.
- `Nc`: control horizon (the control is held constant after `Nc` steps); `move_block`: general
  move-blocking pattern. The exact equivalence requires the defaults (`Nc = Np`, no blocking),
  since blocked moves cannot reproduce the PID trajectory.
- General constraints `lb <= Ax x_p + Au u + Ar r <= ub` on any linear combination of the plant
  state, the control and the reference (`Ax` is m × n_p, `Au` and `Ar` are m × 1), with
  `constraint_ks`, `constraint_soft`, `constraint_binary`, `constraint_prio` and `soft_weight` as
  in `LinearMPCController`. Bounds on the *plant output* are the special case `Ax = C`: the name
  `y` already denotes the controller output (the PID convention of `DiscretePIDStandard`, where
  the limits `y_min`/`y_max` clamp it), so plant-output bounds are written through `Ax` rather
  than through `LinearMPCController`'s `ymin`/`ymax`.
- `binary_controls = [1]` with finite `y_min`/`y_max` restricts the control to the two limit
  values (binary / on-off control, relaxable after `Nc_binary` steps).
- Measured disturbances: `Bd` (n_p × nd) adds a disturbance input to the plant model and the `d`
  input to the component; the model, the constraint predictions and the plant Kalman filter all
  use it, and `disturbance_preview = true` takes the whole known trajectory over the horizon
  (nd*Np column-major entries) instead of the current value. The unconstrained behaviour remains
  *exactly* the PID — the completed-square cost is zero along the PID law whatever the
  disturbance does — so the channel changes the response only where a constraint is active: the
  QP then rejects (and, with preview, pre-acts on) the disturbance where a PID can only react
  through the error. A `u_ff`-style feed-forward that shifts the unconstrained law itself is not
  exactly representable in the QP cost for `Np > 1` and is intentionally not provided.

# Notes
- "Equivalent in the absence of saturation" means the whole *predicted* input trajectory is
  unconstrained: near the limits the MPC departs from the PID as soon as any predicted move would
  saturate, which is what gives it its anti-windup behaviour.
- Without a plant model (`A` left empty) the component reduces to the clamped PID without
  anti-windup — supplying at least a rough plant model is recommended whenever the limits can
  become active.
- All parameters are structural (they are folded into the QP when the model is built) and hence
  not tunable at runtime. `Ts` must equal the period of the clock the component is connected to.
- A PID that does not stabilize the supplied plant model is reported with a warning: the
  unconstrained equivalence still holds, but the constrained behaviour carries no guarantees.
"""
component PIDEquivalentMPCController
  "The reference signal to the controller"
  u_s = RealInput@[clk]() {
    "Dyad": {"placement": {"diagram": {"x1": -100, "y1": 200, "x2": 0, "y2": 300}}}
  }
  "The measurement feedback"
  u_m = RealInput@[clk]() {
    "Dyad": {"placement": {"diagram": {"x1": -100, "y1": 700, "x2": 0, "y2": 800}}}
  }
  "Measured disturbances; nd entries, or nd*Np column-major entries with disturbance preview (only when Bd is given)"
  d = [RealInput@[clk]() for i in 1:nd_in] if has_dist {
    "Dyad": {
      "placement": {"diagram": {"x1": 450, "y1": -100, "x2": 550, "y2": 0, "rot": 90}}
    }
  }
  "The control signal output"
  y = RealOutput@[clk]() {
    "Dyad": {"placement": {"diagram": {"x1": 1000, "y1": 700, "x2": 1100, "y2": 800}}}
  }
  "Solver status of the last step (DAQP exit flag: 1 optimal, 2 with active soft constraints, negative on failure); lmpc_status_name names it"
  exitflag = RealOutput@[clk]() {
    "Dyad": {"placement": {"diagram": {"x1": 1000, "y1": 200, "x2": 1100, "y2": 300}}}
  }
  "Constrained optimizer solving the completed-square QP"
  optimizer = MPCComponents.LinearMPCOptimizer@[clk = clk](mpc = mpc, backend = backend, precision = precision, output_trajectories = output_trajectories, output_statistics = output_statistics) {
    "Dyad": {"placement": {"diagram": {"x1": 400, "y1": 400, "x2": 600, "y2": 600}}}
  }
  "Plant state matrix (n_p × n_p); empty for no plant model (clamped PID without anti-windup)"
  structural parameter A::Real[:, :] = fill(0.0, 0, 0)
  "Plant input matrix (n_p × 1)"
  structural parameter B::Real[:, :] = fill(0.0, 0, 0)
  "Plant output (measurement) matrix (1 × n_p); empty for the identity when n_p = 1"
  structural parameter C::Real[:, :] = fill(0.0, 0, 0)
  "Plant input matrix of the measured disturbances (n_p × nd); empty for none"
  structural parameter Bd::Real[:, :] = fill(0.0, 0, 0)
  "Interpret A and B as continuous-time and ZOH-discretize them with Ts"
  structural parameter continuous::Boolean = false
  "Sample time; must equal the period of the clock the component is connected to"
  structural parameter Ts::Real = 1.0
  "Prediction horizon (only affects the behaviour when the output limits become active)"
  structural parameter Np::Integer = 10
  "Control horizon; the control is held constant after Nc steps. The exact PID equivalence requires the default Nc = Np"
  structural parameter Nc::Integer = Np
  "Take a measured-disturbance trajectory over the horizon (nd*Np entries) instead of the current value"
  structural parameter disturbance_preview::Boolean = false
  "Whether or not to include the integral part"
  structural parameter with_I::Boolean = true
  "Whether or not to include the derivative part"
  structural parameter with_D::Boolean = true
  "Proportional gain"
  structural parameter K::Real = 1.0
  "Integral time constant (only active if `with_I = true`)"
  structural parameter Ti::Real = 1.0
  "Derivative time (only active if `with_D = true`)"
  structural parameter Td::Real = 1.0
  "Maximum derivative gain (only active if `with_D = true`). Typically set to 10-100."
  structural parameter Nd::Real = 10.0
  "Set-point weighting in the proportional part. Set to `0` to prevent step changes in the output due to step changes in the reference."
  structural parameter wp::Real = 1.0
  "Set-point weighting in the derivative part. Set to `0` to prevent very large impulsive changes in the output due to step changes in the reference."
  structural parameter wd::Real = 1.0
  "Anti-windup tracking gain, as in DiscretePIDStandard: the integrator error becomes eI = e + (y - v) / (K Ni), where y is the applied control and v the unconstrained PID law, so the correction is zero whenever no constraint is active (the exact PID equivalence is untouched) and drains the integrator whenever the QP deviates from the law. Inf (the default, unlike DiscretePIDStandard's sqrt(Td/Ti)) disables it: enabling it is recommended whenever output/state constraints are used — without it the integrator winds up against them — but keep it disabled when cost-modifying options (Qy, Qfy, Q3, Nc < Np, move_block) are active, since those deviate from the PID law by design and the correction would fight them. Deliberately kept out of the QP's prediction model, which would let the optimizer steer the integrator and relax its own target law"
  structural parameter Ni::Real = Inf
  "Maximum controller output"
  structural parameter y_max::Real = Inf
  "Minimum controller output"
  structural parameter y_min::Real = ifelse(y_max > 0, -y_max, -Inf)
  "Weight of the completed-square cost; any positive value gives the same unconstrained controller, the value scales the QP"
  structural parameter R_cost::Real = 1.0
  "Stage weight on the squared plant-output tracking error (C x_p - r)^2, added to the PID cost; 0 preserves the exact PID equivalence (requires a plant model)"
  structural parameter Qy::Real = 0.0
  "Terminal weight on the squared plant-output tracking error; 0 preserves the exact PID equivalence (requires a plant model)"
  structural parameter Qfy::Real = 0.0
  "Weight on the squared control increment; 0 preserves the exact PID equivalence"
  structural parameter Q3::Real = 0.0
  "Plant-state coefficient of the general constraints lb <= Ax x_p + Au u + Ar r <= ub (m × n_p); empty for zero. Plant-output bounds are the special case Ax = C"
  structural parameter Ax::Real[:, :] = fill(0.0, 0, 0)
  "Control coefficient of the general constraints (m × 1); empty for zero"
  structural parameter Au::Real[:, :] = fill(0.0, 0, 0)
  "Reference coefficient of the general constraints (m × 1); empty for zero"
  structural parameter Ar::Real[:, :] = fill(0.0, 0, 0)
  "Lower bounds of the general constraints (m); empty for none"
  structural parameter lb::Real[:] = fill(0.0, 0)
  "Upper bounds of the general constraints (m); empty for none"
  structural parameter ub::Real[:] = fill(0.0, 0)
  "Prediction steps on which the general constraints are enforced; empty for 2:Np"
  structural parameter constraint_ks::Integer[:] = []
  "Soften the general constraints"
  structural parameter constraint_soft::Boolean = false
  "Enforce the general constraints as binary (one of the bounds holds with equality)"
  structural parameter constraint_binary::Boolean = false
  "Priority of the general constraints"
  structural parameter constraint_prio::Integer = 0
  "Penalty weight on soft-constraint violations"
  structural parameter soft_weight::Real = 1e6
  "DAQP primal feasibility tolerance; -1 keeps the solver default. Applied to the Julia solver and embedded in the generated C"
  structural parameter primal_tol::Real = -1.0
  "DAQP dual feasibility (optimality) tolerance; -1 keeps the solver default"
  structural parameter dual_tol::Real = -1.0
  "DAQP iteration limit; -1 keeps the solver default"
  structural parameter iter_limit::Integer = -1
  "Indices of the controls restricted to the values {y_min, y_max}; empty for none. The only valid index is 1"
  structural parameter binary_controls::Integer[:] = []
  "Number of initial horizon steps over which the binary controls are restricted to {y_min, y_max}, the rest are relaxed to continuous; -1 for the whole control horizon"
  structural parameter Nc_binary::Integer = -1
  "Move-blocking pattern, the number of steps each control move is held; empty for none. Non-default blocking breaks the exact PID equivalence"
  structural parameter move_block::Integer[:] = []
  "Process-noise covariance of the plant Kalman filter (n_p × n_p, or n_p diagonal entries); empty for the identity"
  structural parameter R1::Real[:, :] = fill(0.0, 0, 0)
  "Measurement-noise covariance of the plant Kalman filter (1 × 1); empty for the identity"
  structural parameter R2::Real[:, :] = fill(0.0, 0, 0)
  "QP solver backend; Auto uses the compiled C code when a C compiler is available"
  structural parameter backend::Backend = Backend.Auto()
  "Float type of the generated QP code (C backend only)"
  structural parameter precision::MPCPrecision = MPCPrecision.Double()
  "Also record the predicted state and control trajectories and the QP statistics of every tick (variables optimizer.x_pred, optimizer.u_pred, optimizer.stats), for debugging and mpc_gui; Julia backend only. The readbacks are three equations whose number does not grow with the prediction horizon, though the trajectories themselves are stored in the solution; output_statistics records the statistics alone"
  structural parameter output_trajectories::Boolean = false
  "Record the QP statistics of every tick (variable optimizer.stats) without the predicted trajectories; Julia backend only. One equation. Implied by output_trajectories"
  structural parameter output_statistics::Boolean = false
  "Plant state dimension (derived)"
  final structural parameter np::Integer = size(A, 1)
  "Number of measured disturbances (derived)"
  final structural parameter nd::Integer = size(Bd, 2)
  "Whether measured disturbances are present (derived)"
  final structural parameter has_dist::Boolean = nd > 0
  "Width of the disturbance input: nd, or nd*Np with disturbance preview (derived)"
  final structural parameter nd_in::Integer = ifelse(disturbance_preview, nd * Np, nd)
  "Width of the internal disturbance vector, max(nd_in, 1) (derived)"
  final structural parameter nd_in_eff::Integer = max(nd_in, 1)
  "Whether a plant model was given (derived)"
  final structural parameter has_plant::Boolean = np > 0
  "Width of the plant-estimator vectors, max(np, 1) (derived)"
  final structural parameter np_eff::Integer = max(np, 1)
  "State dimension of the augmented model [x_p; I; D; r; y] (derived)"
  final structural parameter nxt::Integer = np + ifelse(with_I, 1, 0) + ifelse(with_D, 1, 0) + 2
  "Index of the integrator state I in the augmented state (derived)"
  final structural parameter i_I::Integer = np + 1
  "Index of the derivative-filter state D in the augmented state (derived)"
  final structural parameter i_D::Integer = np + ifelse(with_I, 1, 0) + 1
  "Index of the constant-reference state r in the augmented state (derived)"
  final structural parameter i_r::Integer = nxt - 1
  "Index of the measured-output state y in the augmented state (derived)"
  final structural parameter i_y::Integer = nxt
  "Discrete-time plant state matrix, padded to np_eff × np_eff (derived)"
  final structural parameter F::Real[np_eff, np_eff] = _pidmpc_F(A, B, C, continuous, Ts)
  "Discrete-time plant input matrix, padded to np_eff × 1 (derived)"
  final structural parameter G::Real[np_eff, 1] = _pidmpc_G(A, B, C, continuous, Ts)
  "Plant output matrix, padded to 1 × np_eff (derived)"
  final structural parameter Cp::Real[1, np_eff] = _pidmpc_C(A, B, C, continuous, Ts)
  "Discrete-time plant disturbance input matrix, padded to np_eff × nd_in_eff; only the first nd columns (the current disturbance) are nonzero (derived)"
  final structural parameter Gd::Real[np_eff, nd_in_eff] = _pidmpc_Gd(A, B, C, Bd, continuous, Ts, nd_in_eff)
  "Augmented model state matrix (derived)"
  final structural parameter Atilde::Real[nxt, nxt] = _pidmpc_A(F, G, Cp, np, Ts, Ti, Td, Nd, with_I, with_D)
  "Augmented model input matrix (derived)"
  final structural parameter Btilde::Real[nxt, 1] = _pidmpc_B(G, Cp, np, Ts, Td, Nd, with_I, with_D)
  "Augmented model disturbance input matrix (derived)"
  final structural parameter Btilde_d::Real[nxt, nd] = _pidmpc_Bd(Gd, Cp, np, nd, Ts, Td, Nd, with_I, with_D)
  "The PID law as state feedback u = -Ltilde x, Ltilde = [0, -K, -K, -K wp, K] (derived)"
  final structural parameter Ltilde::Real[1, nxt] = _pidmpc_L(Atilde, Btilde, np, K, wp, with_I, with_D)
  "Row extracting the plant-output tracking error C x_p - r from the augmented state (derived)"
  final structural parameter Ttilde::Real[1, nxt] = _pidmpc_tracking_row(Cp, np, with_I, with_D, Qy, Qfy)
  "Stage state weight of the QP: the completed square R_cost (u + Ltilde x)^2 contributes R_cost Ltilde' Ltilde, the optional plant-output tracking cost Qy (C x_p - r)^2 adds Qy Ttilde' Ttilde (derived)"
  final structural parameter Q1t::Real[nxt, nxt] = R_cost * transpose(Ltilde) * Ltilde + Qy * transpose(Ttilde) * Ttilde
  "State-control cross weight of the completed square, R_cost Ltilde' (derived)"
  final structural parameter Q12t::Real[nxt, 1] = R_cost * transpose(Ltilde)
  "Control weight of the completed square, R_cost (derived)"
  final structural parameter Q2t::Real[1, 1] = fill(R_cost, 1, 1)
  "Terminal state weight, the plant-output tracking cost Qfy (C x_p - r)^2 at the end of the horizon (derived)"
  final structural parameter Qfxt::Real[nxt, nxt] = Qfy * transpose(Ttilde) * Ttilde
  "LinearMPC.jl controller with the completed-square cost (derived)"
  final structural parameter mpc::MPCRef = build_linear_mpc(A = Atilde, B = Btilde, Bd = Btilde_d, disturbance_preview = disturbance_preview, Ts = Ts, Np = Np, Nc = Nc, Q1 = Q1t, Q2 = Q2t, Q3 = _pidmpc_Q3(Q3), Q12 = Q12t, Qf = _pidmpc_Qf(nxt), Qfx = Qfxt, reference_tracking = false, umin = _pidmpc_umin(y_min, y_max), umax = _pidmpc_umax(y_min, y_max), Ax = _pidmpc_Ax(Ax, Ar, np, with_I, with_D), Au = Au, lb = lb, ub = ub, constraint_ks = constraint_ks, constraint_soft = constraint_soft, constraint_binary = constraint_binary, constraint_prio = constraint_prio, soft_weight = soft_weight, primal_tol = primal_tol, dual_tol = dual_tol, iter_limit = iter_limit, binary_controls = binary_controls, Nc_binary = Nc_binary, move_block = move_block, state_labels = _pidmpc_state_labels(np, with_I, with_D), control_labels = ["y"])
  "Initial plant state estimate"
  parameter x0::Real[np_eff] = fill(0.0, np_eff)
  "Derivative-filter pole, as in the backward-difference filter of DiscretePIDStandard (derived)"
  final parameter a_d::Real = Td / (Td + Nd * Ts)
  "Derivative-filter gain (derived)"
  final parameter b_d::Real = Nd * Td / (Td + Nd * Ts)
  "Plant estimator state matrix, F as a runtime parameter (derived)"
  final parameter F_p::Real[np_eff, np_eff] = F
  "Plant estimator input matrix, G as a runtime parameter (derived)"
  final parameter G_p::Real[np_eff, 1] = G
  "Plant estimator measurement matrix, Cp as a runtime parameter (derived)"
  final parameter C_p::Real[1, np_eff] = Cp
  "Plant estimator disturbance input matrix, Gd as a runtime parameter (derived)"
  final parameter Gd_p::Real[np_eff, nd_in_eff] = Gd
  "Steady-state Kalman gain of the plant estimator (derived)"
  final parameter K_kf::Real[np_eff, 1] = _pidmpc_kf_K(F, G, Cp, np, R1, R2)
  "Zero placeholder for the disturbance vector when d is absent"
  final parameter d_dummy::Real[nd_in_eff] = fill(0.0, nd_in_eff)
  "Error signal"
  variable e::Real
  "Setpoint-weighted error for the derivative"
  variable wde::Real
  "State of the integrator"
  variable I_state::Real
  "State of the filtered derivative"
  variable D_state::Real
  "Corrected plant state estimate"
  variable xcorr::Real[np_eff]
  "Predicted plant state estimate"
  variable xpred::Real[np_eff]
  "Measurement as a vector, for the estimator equations"
  variable ym::Real[1]
  "Applied control as a vector, for the estimator equations"
  variable uv::Real[1]
  "Measured disturbances as a vector, for the estimator equations"
  variable d_int::Real[nd_in_eff]
  "Unconstrained PID law, for the anti-windup tracking correction"
  variable v_pid::Real
  "Integrator error including the anti-windup tracking signal"
  variable eI::Real
relations
  e = u_s - u_m
  wde = wd * u_s - u_m
  v_pid = K * (wp * u_s - u_m + I_state + D_state)
  eI = e + (y - v_pid) / (K * Ni)
  ym[1] = u_m
  uv[1] = y
  # The augmented state [x_p; I; D; r; y] handed to the optimizer
  if with_I
    initial I_state@(clk-1) = 0
    initial eI@(clk-1) = 0
    I_state@clk = I_state@(clk-1) + Ts / Ti * eI@(clk-1)
    optimizer.xhat[i_I] = I_state
  else
    I_state = 0
  end
  if with_D
    initial D_state@(clk-1) = 0
    initial wde@(clk-1) = 0
    D_state@clk = a_d * D_state@(clk-1) + b_d * (wde@clk - wde@(clk-1))
    optimizer.xhat[i_D] = D_state
  else
    D_state = 0
  end
  optimizer.xhat[i_r] = u_s
  optimizer.xhat[i_y] = u_m
  if has_dist
    d_int = d
  else
    d_int = d_dummy
  end
  if has_plant
    initial xpred@(clk-1) = x0
    xcorr@clk = xpred@(clk-1) + K_kf * (ym@clk - C_p * xpred@(clk-1))
    xpred@clk = F_p * xcorr@clk + G_p * uv@clk + Gd_p * d_int@clk
    for i in 1:np
      optimizer.xhat[i] = xcorr[i]
    end
  else
    xcorr = fill(0.0, np_eff)
    xpred = fill(0.0, np_eff)
  end
  if has_dist
    for i in 1:nd_in
      connect(d[i], optimizer.d[i])
    end
  end
  connect(optimizer.u[1], y) {}
  connect(optimizer.exitflag, exitflag) {}
metadata {
  "Dyad": {
    "icons": {"default": "dyad://MPCComponents/PIDEquivalentMPCController.svg"},
    "doc": {"behavior": false},
    "panel": {
      "control": "tabs",
      "contents": [
        "K",
        "Ti",
        "Td",
        "Nd",
        "wp",
        "wd",
        "Ni",
        "y_max",
        "y_min",
        "with_I",
        "with_D",
        "Ts"
      ],
      "children": {
        "Plant Model": {"contents": ["A", "B", "C", "Bd", "continuous"]},
        "State Estimator": {"contents": ["R1", "R2", "x0"]},
        "MPC": {
          "children": {
            "Horizon": {"contents": ["Np", "Nc", "move_block", "disturbance_preview"]},
            "Cost": {"contents": ["R_cost", "Qy", "Qfy", "Q3"]},
            "Constraints": {
              "contents": [
                "Ax",
                "Au",
                "Ar",
                "lb",
                "ub",
                "constraint_ks",
                "constraint_soft",
                "constraint_binary",
                "constraint_prio",
                "soft_weight",
                "binary_controls",
                "Nc_binary"
              ]
            }
          }
        },
        "Solver": {
          "contents": [
            "backend",
            "precision",
            "output_trajectories",
            "output_statistics",
            "primal_tol",
            "dual_tol",
            "iter_limit"
          ]
        },
        "Internal": {
          "contents": [
            "np",
            "nd",
            "has_dist",
            "nd_in",
            "nd_in_eff",
            "has_plant",
            "np_eff",
            "nxt",
            "i_I",
            "i_D",
            "i_r",
            "i_y",
            "F",
            "G",
            "Cp",
            "Gd",
            "Btilde_d",
            "d_dummy",
            "d_int",
            "Atilde",
            "Btilde",
            "Ltilde",
            "Ttilde",
            "Q1t",
            "Q12t",
            "Q2t",
            "Qfxt",
            "mpc",
            "a_d",
            "b_d",
            "v_pid",
            "eI",
            "F_p",
            "G_p",
            "C_p",
            "Gd_p",
            "K_kf",
            "e",
            "wde",
            "I_state",
            "D_state",
            "xcorr",
            "xpred",
            "ym",
            "uv"
          ]
        }
      }
    }
  }
}
end

Test Cases

No test cases defined.

Related

• Examples

• Experiments

• Analyses