ANALYSIS

PID Autotuning Analysis

Download as a Dyad projectPIDAnalysis.zipOpen in Dyad Studio

The PID autotuning analysis is used to automatically tune the parameters of a PID controller for a given model. The analysis computes optimal PID gains that satisfy frequency-domain robustness constraints, such as maximum sensitivity (Ms), complementary sensitivity (Mt), peak noise sensitivity (Mks), and integrated noise sensitivity (Mks2).

Method Overview

PID autotuning in DyadControlSystems is based on step-response optimization, subject to robustness constraints on closed-loop sensitivity functions. The method performs joint optimization of a PID controller and a measurement filter on either of the forms

by solving

subject to

where is the control error and is a user-defined metric (defaults to abs2). The user can select which signals to optimize, the reference value, and the filter order.

The transfer functions in the frequency-domain constraints are defined as follows:

Sensitivity:
Complementary sensitivity:
Noise sensitivity:

See control-system sensitivity and DyadControlSystems autotuning documentation for more information.

Example Definition

The following example will demonstrate autotuning of a PID controller for an active suspension system. The model of the system is available in DyadExampleComponents.

Since the PID autotuning analysis internally opens the loop at the measurement point y, the connection from seat_pos.s to pid.u_m is broken during linearization, and pid.u_m becomes a free input that needs an operating-point value. We also need to pin the algebraic spring-damper forces and the internal PID-block outputs so that the linearization solver does not encounter cyclic guesses. We use Dyad's guess keyword to provide these operating-point values:

julia

using DyadExampleComponents, DyadInterface
ActiveSuspensionPIDAutotuningAnalysis() # hide
solution = ActiveSuspensionPIDAutotuningAnalysis()

┌ Warning: The initialization system is structurally singular. Guess values may significantly affect the initial values of the ODE. The problematic variables are SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}[car_and_suspension₊mass₊s(t), car_and_suspension₊spring_damper₊s_rel0, seat₊spring_damper₊s_rel0].
│
│ Note that the identification of problematic variables is a best-effort heuristic.
└ @ ModelingToolkitBase ~/.julia/packages/ModelingToolkitBase/N5KVH/src/problems/initializationproblem.jl:140
┌ Warning: Initialization system is underdetermined. 0 equations for 11 unknowns. Initialization will default to using least squares. `SCCNonlinearProblem` can only be used for initialization of fully determined systems and hence will not be used here.
│
│
│ Call `analyze_initialization_jacobian(prob)` on the constructed problem to see which unknowns are underdetermined (and which equations, if any, are redundant).
│
│ To suppress this warning, pass `warn_initialize_determined = false`. To turn this warning into an error, pass `fully_determined = true`.
└ @ ModelingToolkit ~/.julia/packages/ModelingToolkit/ZgEKM/src/initialization.jl:25

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit https://github.com/coin-or/Ipopt
******************************************************************************

This is Ipopt version 3.14.19, running with linear solver MUMPS 5.9.0.

Number of nonzeros in equality constraint Jacobian...:        0
Number of nonzeros in inequality constraint Jacobian.:     2404
Number of nonzeros in Lagrangian Hessian.............:        0

Total number of variables............................:        4
                     variables with only lower bounds:        4
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:        0
Total number of inequality constraints...............:      601
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:      601

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  1.1089185e-02 1.45e+01 1.14e-01   0.0 0.00e+00    -  0.00e+00 0.00e+00   0
   5  1.0522811e-02 1.29e+01 2.27e-01  -5.4 6.70e+01    -  6.31e-03 3.31e-03h  1
  10  7.4366640e-03 2.28e+03 2.13e+01  -0.7 1.78e+03    -  2.15e-02 1.27e-02f  1
  15  7.1678939e-03 0.00e+00 1.20e+02   0.2 3.77e+02    -  2.38e-02 2.44e-01f  1
  20  9.2469686e-03 0.00e+00 5.67e+01  -0.0 5.18e+01    -  8.86e-01 1.00e+00h  1
  25  5.7540576e-03 0.00e+00 7.86e-03  -4.2 1.60e+01    -  9.92e-01 1.00e+00h  1
  30  4.1134438e-03 0.00e+00 1.22e-02  -4.5 5.94e+01    -  1.00e+00 1.00e+00h  1
  35  1.0104852e-03 0.00e+00 5.79e-04  -5.4 2.74e+01    -  9.99e-01 1.00e+00h  1
  40  6.0712668e-04 0.00e+00 7.54e-04  -6.7 2.39e+00    -  6.91e-01 1.00e+00h  1

Number of Iterations....: 42

                                   (scaled)                 (unscaled)
Objective...............:   6.0118943309683435e-04    6.0118943309683435e-04
Dual infeasibility......:   2.0811264791115461e-07    2.0811264791115461e-07
Constraint violation....:   0.0000000000000000e+00    0.0000000000000000e+00
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   1.1377472139212053e-08    1.1377472139212053e-08
Overall NLP error.......:   2.0811264791115461e-07    2.0811264791115461e-07


Number of objective function evaluations             = 43
Number of objective gradient evaluations             = 43
Number of equality constraint evaluations            = 0
Number of inequality constraint evaluations          = 44
Number of equality constraint Jacobian evaluations   = 0
Number of inequality constraint Jacobian evaluations = 43
Number of Lagrangian Hessian evaluations             = 0
Total seconds in IPOPT                               = 103.203

EXIT: Optimal Solution Found.
┌ Warning: The initialization system is structurally singular. Guess values may significantly affect the initial values of the ODE. The problematic variables are SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}[car_and_suspension₊mass₊s(t), car_and_suspension₊spring_damper₊s_rel0, seat₊spring_damper₊s_rel0].
│
│ Note that the identification of problematic variables is a best-effort heuristic.
└ @ ModelingToolkitBase ~/.julia/packages/ModelingToolkitBase/N5KVH/src/problems/initializationproblem.jl:140
┌ Warning: Initialization system is underdetermined. 0 equations for 11 unknowns. Initialization will default to using least squares. `SCCNonlinearProblem` can only be used for initialization of fully determined systems and hence will not be used here.
│
│
│ Call `analyze_initialization_jacobian(prob)` on the constructed problem to see which unknowns are underdetermined (and which equations, if any, are redundant).
│
│ To suppress this warning, pass `warn_initialize_determined = false`. To turn this warning into an error, pass `fully_determined = true`.
└ @ ModelingToolkit ~/.julia/packages/ModelingToolkit/ZgEKM/src/initialization.jl:25
This is Ipopt version 3.14.19, running with linear solver MUMPS 5.9.0.

Number of nonzeros in equality constraint Jacobian...:        0
Number of nonzeros in inequality constraint Jacobian.:     2404
Number of nonzeros in Lagrangian Hessian.............:        0

Total number of variables............................:        4
                     variables with only lower bounds:        4
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:        0
Total number of inequality constraints...............:      601
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:      601

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  1.1089185e-02 1.45e+01 1.14e-01   0.0 0.00e+00    -  0.00e+00 0.00e+00   0
   5  1.0522811e-02 1.29e+01 2.27e-01  -5.4 6.70e+01    -  6.31e-03 3.31e-03h  1
  10  7.4366640e-03 2.28e+03 2.13e+01  -0.7 1.78e+03    -  2.15e-02 1.27e-02f  1
  15  7.1678939e-03 0.00e+00 1.20e+02   0.2 3.77e+02    -  2.38e-02 2.44e-01f  1
  20  9.2469686e-03 0.00e+00 5.67e+01  -0.0 5.18e+01    -  8.86e-01 1.00e+00h  1
  25  5.7540576e-03 0.00e+00 7.86e-03  -4.2 1.60e+01    -  9.92e-01 1.00e+00h  1
  30  4.1134438e-03 0.00e+00 1.22e-02  -4.5 5.94e+01    -  1.00e+00 1.00e+00h  1
  35  1.0104852e-03 0.00e+00 5.79e-04  -5.4 2.74e+01    -  9.99e-01 1.00e+00h  1
  40  6.0712668e-04 0.00e+00 7.54e-04  -6.7 2.39e+00    -  6.91e-01 1.00e+00h  1

Number of Iterations....: 42

                                   (scaled)                 (unscaled)
Objective...............:   6.0118943309683435e-04    6.0118943309683435e-04
Dual infeasibility......:   2.0811264791115461e-07    2.0811264791115461e-07
Constraint violation....:   0.0000000000000000e+00    0.0000000000000000e+00
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   1.1377472139212053e-08    1.1377472139212053e-08
Overall NLP error.......:   2.0811264791115461e-07    2.0811264791115461e-07


Number of objective function evaluations             = 43
Number of objective gradient evaluations             = 43
Number of equality constraint evaluations            = 0
Number of inequality constraint evaluations          = 44
Number of equality constraint Jacobian evaluations   = 0
Number of inequality constraint Jacobian evaluations = 43
Number of Lagrangian Hessian evaluations             = 0
Total seconds in IPOPT                               = 1.485

EXIT: Optimal Solution Found.

julia

using Plots
figs  = artifacts(solution, :SensitivityFunctions)
fign  = artifacts(solution, :NoiseSensitivityAndController)
figr  = artifacts(solution, :OptimizedResponse)
figny = artifacts(solution, :NyquistPlot)

plot!(figs, ylims=(1e-2, Inf))
plot!(figny, legend=:bottomleft)

plot(figs, fign, figr, figny,
    link = :none,
    size = (600, 600),
)

We may also ask for the optimized controller parameters:

julia

artifacts(solution, :OptimizedParameters)

1×8 DataFrame

Row	kp_parallel	ki_parallel	kd_parallel	Kp_standard	Ti_standard	Td_standard	Tf	Nd
	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64
1	1358.06	7573.58	854.719	1358.06	0.179316	0.629365	0.0298335	21.0959

Using the Tuned Controller in a Dyad Model

Use apply to load the tuned parameters from the :ControllerParameters TOML artifact:

julia

print(artifacts(solution, :ControllerParameters))

# PID controller parameters tuned by PIDAutotuningAnalysis.
# Place this file in the directory `assets` of a Dyad package `MyPackage` and load it with `apply`, e.g.,
#   controller = BlockComponents.Continuous.LimPID(apply "dyad://MyPackage/pid.toml")
# Structural parameters set explicitly in the modification take precedence over `with_I` and `with_D` in this file.

with_I = true
with_D = true
k = 1358.0646941911252
Ti = 0.17931613230860868
Td = 0.6293654099068342
Nd = 21.095924191625503

To save the parameters, add output_file = "pid.toml" to the extends DyadControlSystems.PIDAutotuningAnalysis(...) argument list in your analysis definition. Running the analysis writes pid.toml to the working directory. Move it to assets/pid.toml in MyPackage, then load it:

dyad

controller = BlockComponents.Continuous.LimPID(apply "dyad://MyPackage/pid.toml")

The file's header comments identify the required controller structure:

Analysis settings	Model structure
Continuous, `filter_order = 1` (default)	`LimPID` realizes the tuned controller exactly.
Continuous, `filter_order = 2`	Add `BlockComponents.Continuous.SecondOrder` low-pass filters named `filter_y` and `filter_r` for the measurement and reference. Apply the file to their parent component, which must also contain a `LimPID` named `pid`. The file has a table for each subcomponent.
`discrete_time = true`	Use `DiscreteComponents.DiscretePIDStandard` with sample interval `Ts`.

Tuning Several Controllers

PID autotuning handles one controller in one feedback loop. For broader tuning problems, use DyadControlSystems' StructuredAutoTuningProblem from Julia. It supports:

Multiple controllers, including cascade loops and interacting PIDs.
Other controller structures.
Multiple operating points or uncertain parameters.

It tunes arbitrary model parameters against time- and frequency-domain objectives at analysis points. It is not available as a Dyad analysis.

Export the result with parameter_file(res, path). The TOML file uses model-relative parameter names, with one table per subcomponent:

toml

[inner_controller]
k = 2.1
Ti = 0.12

[outer_controller]
k = 8.5

Apply it to the component used to create the tuned model:

dyad

model = MyModel(apply "dyad://MyPackage/parameters.toml")

Analysis Arguments

The following arguments define a PID autotuning analysis:

Required Arguments

model: The system model to be controlled.
measurement::String: Name of the measured output signal.
control_input::String: Name of the control input signal.
step_input::String: Name of the step input signal.
step_output::String: Name of the step output signal.
Ts::Real: Discretization sampling time; also the controller sample interval when discrete_time = true.
duration::Real: Simulation duration.
Ms::Real: Maximum sensitivity constraint. Values between 1 and 2 are typical. A lower value increases robustness but may reduce performance.
Mt::Real: Maximum complementary sensitivity constraint. Values between 1 and 2 are typical. A lower value increases robustness but may reduce performance.
Mks::Real: Maximum noise sensitivity constraint. This controls the maximum amplification of measurement noise in the control signal.

Optional Arguments

ref::Real = 0.0: Reference value for the step input.
Mks2::Real = 1e300: Upper bound on the H2 norm of the transfer function from measurement noise to the control signal. The default disables this constraint. It can be used independently of Mks or together with it.
exact_h2::Bool = false: Whether to evaluate the Mks2 constraint using a Lyapunov equation over the full frequency range for continuous-time tuning. With false, the analysis integrates over the frequency grid. See Integrated noise sensitivity.
disc::String = "tustin": Discretization method.
kp_lb, ki_lb, kd_lb, Tf_lb: Lower bounds for PID parameters.
kp_ub, ki_ub, kd_ub, Tf_ub: Upper bounds for PID parameters.
kp_guess, ki_guess, kd_guess, Tf_guess: Initial guesses for PID parameters.
timeweight::Bool = false: Whether or not to use time-weighted optimization.
filter_order::Int = 1: Noise filter order (1 or 2). 1 filters only the derivative term; 2 filters the measurement and reference, acting on the entire PID controller.
optimize_d::Bool = false: Whether or not to optimize the damping parameter of the second-order filter.
discrete_time::Bool = false: Tune a DiscreteComponents.DiscretePIDStandard with sample interval Ts and default discretization methods. Requires filter_order = 1.
output_file::String = "": Path for the tuned controller TOML file. Empty disables writing.
wl::Real: Lower frequency bound for design.
wu::Real: Upper frequency bound for design.
num_frequencies::Int: Number of frequency points.
soft::Bool: Whether or not to use soft constraints (may be beneficial for difficult problems).
soft_penalty::Float64: The penalty factor for violating soft constraints. Defaults to 1e6.
exact_hessian::Bool = false: Whether or not to compute the exact Hessian using AD. This may improve performance in some cases, in particular if the solver requires many steps to converge. If exact Hessian is turned off, a Hessian-free algorithm will be used (IPOPT BFGS approximation).

Controller Structure

The optimized controller has the form:

for filter_order = 2, where is a second-order low-pass filter. With filter_order = 1, only the derivative term is filtered, as shown in Method Overview.

Constraints

Sensitivity (Ms): Limits the peak of the sensitivity function . Typical values are between 1 and 2.
Complementary Sensitivity (Mt): Limits the peak of . Typical values are between 1 and 2.
Noise Sensitivity (Mks): Limits the peak of , which determines the amplification of measurement noise in the control signal.
Integrated Noise Sensitivity (Mks2): Limits the H2 norm of , which accounts for the amplification of measurement noise across frequencies.

Lower values increase robustness but may reduce performance.

The peak in the sensitivity function is associated with lower bounds on the classical gain and phase margins, :

For example, a maximum sensitivity of corresponds to a gain margin of at least 4.3 and a phase margin of at least 45 degrees. The peak constraints and corresponds to circles which the Nyquist curve must fall outside, these circles are drawn with dashed lines in the Nyquist plot artifact.

Integrated noise sensitivity

The peak constraint Mks permits a nearly constant noise gain above the filter cutoff frequency. The H2 constraint Mks2 limits the integral of the squared gain, so a controller with high noise gain over a broad frequency interval requires a larger bound. For a stable continuous-time closed loop,

For unit-intensity white measurement noise, this norm determines the standard deviation of the control signal. Mks and Mks2 have different units and numerical values; their bounds are selected independently.

With exact_h2 = false, the analysis applies the trapezoidal rule on the grid determined by wl, wu, and num_frequencies. It integrates with respect to angular frequency, including when the grid is logarithmic. Frequencies outside this interval are excluded. Use a sufficiently broad and dense grid and check the result after increasing the frequency range and resolution. For a first-order derivative filter, a finite integration interval represents a finite measurement-noise bandwidth.

With exact_h2 = true, the analysis evaluates the full-band norm using a Lyapunov equation. This method supports continuous-time tuning and requires filter_order = 2 and a strictly positive Tf_lb. A first-order derivative filter generally leaves a nonzero direct feedthrough from measurement noise to the control signal, which gives an infinite continuous-time H2 norm. A second-order filter acts on the entire PID controller and provides the required high-frequency attenuation. In the Julia interface, a custom lower bound on the damping parameter must also be strictly positive when optimize_d = true.

The following configuration constrains integrated noise sensitivity for the active suspension model:

dyad

analysis ActiveSuspensionPIDH2Analysis
  extends DyadControlSystems.PIDAutotuningAnalysis(
    measurement = "y",
    control_input = "u",
    duration = 10.0,
    Ts = 0.01,
    Ms = 1.2,
    Mt = 1.2,
    Mks2 = 1e5,
    exact_h2 = true,
    filter_order = 2,
    Tf_lb = 1e-4
  )
  model = ActiveSuspensionForAutotuning()
end

The default Mks leaves the peak noise gain unconstrained in this example. Set a finite Mks to impose both bounds. The Julia interface uses the same Mks2 and Boolean exact_h2 options in AutoTuningProblem.

For discrete-time tuning (discrete_time = true), use exact_h2 = false. The frequency approximation uses

The factor matches the usual discrete-time H2 convention for unit-variance white measurement noise.

Advanced Options

Reference vs. Disturbance Step: By default, the disturbance step response is optimized. To optimize reference tracking, set step_input = "reference_input".
Controller Type: By setting parameter bounds, you can restrict the controller to P, PI, PD, or PID.
Filter Damping: With filter_order = 2 and optimize_d = true, the damping ratio of the filter is also optimized.
Penalizing Control Action: You can augment the model to penalize control effort or noise amplification.

Artifacts

A PIDAutotuningAnalysis returns the following artifacts:

Standard Plots

:SensitivityFunctions: Bode plot of sensitivity (S) and complementary sensitivity (T) functions.
:NoiseSensitivityAndController: Bode plot of noise sensitivity (KS) and controller (K).
:OptimizedResponse: Step response of the optimized closed-loop system.
:NyquistPlot: Nyquist plot of the loop transfer function.

Tables

OptimizedParameters: DataFrame of the optimized PID parameters. The parameters are available for the parallel controller structure as well as converted to the standard form

Files

ControllerParameters: Tuned controller parameters as TOML for apply.

Scalar Types

Partial Types

Enumerated Types

Function Types

Partial Types

Enumerated Types

Partial Types

Enumerated Types

Partial Types

Enumerated Types

Partial Types

Scalar Types

Partial Types

Partial Types

Types

Partial Types

Enumerated Types

Struct Types

Connectors

Components

Analyses

Examples

Tests

Scalar Types

Partial Types

Enumerated Types

Enumerated Types

Types

Scalar Types

Partial Types

Enumerated Types

Components

Analyses

Examples

PID Autotuning Analysis ​

Method Overview ​

Example Definition ​

Using the Tuned Controller in a Dyad Model ​

Tuning Several Controllers ​

Analysis Arguments ​

Required Arguments ​

Optional Arguments ​

Controller Structure ​

Constraints ​

Integrated noise sensitivity ​

Advanced Options ​

Artifacts ​

Standard Plots ​

Tables ​

Files ​

Further Reading ​

PID Autotuning Analysis

Method Overview

Example Definition

Using the Tuned Controller in a Dyad Model

Tuning Several Controllers

Analysis Arguments

Required Arguments

Optional Arguments

Controller Structure

Constraints

Integrated noise sensitivity

Advanced Options

Artifacts

Standard Plots

Tables

Files

Further Reading