Chapter 29: Industrial and Mechatronic Case Studies in Linear Control

Lesson 5: End-to-End Linear Controller Design Mini-Project

In this lesson, we carry out an end-to-end SISO linear controller design for a mechatronic actuator (a DC-motor-driven rotary stage), starting from specifications, passing through plant modeling, controller synthesis (time- and frequency-domain), and finishing with multi-language implementations in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica. The goal is to consolidate all classical linear control concepts developed in previous chapters into a single coherent mini-project.

1. Mini-Project Overview and Architecture

We consider a DC-motor-driven rotary stage whose angular position \( \theta(t) \) must track a step reference with tight requirements on transient behavior and steady-state accuracy. The control architecture is a unity-feedback loop with a tunable PD controller \( C(s) \) acting on a plant transfer function \( G(s) \).

Design steps:

  • Specify time-domain and frequency-domain performance requirements.
  • Derive or recall the linear plant model from the DC motor equations.
  • Choose a controller structure (here: PD) consistent with the plant type.
  • Compute controller gains analytically from pole-placement formulas.
  • Verify closed-loop performance via time and frequency response.
  • Implement the controller and simulation in multiple programming languages.
flowchart TD
  S["Specs: ts, Mp, ess, bandwidth"] --> P["Linear model G(s) from DC motor"]
  P --> A["Open-loop analysis: type, poles, zeros"]
  A --> CSTR["Select controller C(s) (here PD)"]
  CSTR --> D["Solve for gains from desired poles"]
  D --> V["Simulate time & frequency response"]
  V --> IMPL["Implement in Python / C++ / Java / MATLAB / Mathematica"]
        

2. Plant Model for the Rotary Stage

The rotary stage is driven by a DC motor with rotor inertia \( J \), viscous friction coefficient \( b \), torque constant \( K_t \), back-EMF constant \( K_e \), armature resistance \( R \), and inductance \( L \). The mechanical and electrical dynamics are

\[ J \ddot{\theta}(t) + b \dot{\theta}(t) = K_t i(t), \]

\[ L \dot{i}(t) + R i(t) + K_e \dot{\theta}(t) = u(t), \]

where \( u(t) \) is the armature voltage, \( i(t) \) is the armature current, and \( \theta(t) \) is the load angle.

Taking Laplace transforms with zero initial conditions yields

\[ J s^2 \Theta(s) + b s \Theta(s) = K_t I(s), \]

\[ L s I(s) + R I(s) + K_e s \Theta(s) = U(s). \]

In many low-power motors, the electrical dynamics are much faster than the mechanical dynamics. Approximating \( L \approx 0 \) (quasi-static current loop) gives

\[ R I(s) + K_e s \Theta(s) = U(s) \;\Rightarrow\; I(s) = \frac{U(s) - K_e s \Theta(s)}{R}. \]

Substituting into the mechanical equation:

\[ J s^2 \Theta(s) + b s \Theta(s) = K_t \frac{U(s) - K_e s \Theta(s)}{R}. \]

Rearranging:

\[ \left(J s^2 + \left(b + \frac{K_t K_e}{R}\right) s \right)\Theta(s) = \frac{K_t}{R} U(s). \]

Hence the transfer function from voltage to angle is

\[ \frac{\Theta(s)}{U(s)} = \frac{\dfrac{K_t}{R}}{s\left(J s + b + \dfrac{K_t K_e}{R}\right)} = \frac{K}{s(T s + 1)}, \]

where we define the equivalent gain and time constant

\[ K = \frac{K_t}{R \left(b + \dfrac{K_t K_e}{R}\right)}, \quad T = \frac{J}{b + \dfrac{K_t K_e}{R}}. \]

For the mini-project, we assume that previous experimental identification (or detailed modeling in Chapter 29 Lesson 1) has given the normalized plant

\[ G(s) = \frac{\Theta(s)}{U(s)} = \frac{1}{s(0.5 s + 1)}. \]

That is, \( K = 1 \), \( T = 0.5 \, \text{s} \). The plant has a pole at the origin and a real pole at \( s = -2 \), and is a type-1 system.

3. Performance Specifications and Desired Closed-Loop Poles

We specify classical time-domain requirements for the closed-loop position response to a unit step reference:

  • Settling time (2% criterion) \( t_s \) with \( t_s \leq 1.0 \, \text{s} \).
  • Overshoot \( M_p \) with \( M_p \leq 10\% \).
  • Zero steady-state error to step input: \( e_{ss} = 0 \).

For a standard second-order system with transfer function

\[ T_{2}(s) = \frac{\omega_n^2}{s^2 + 2\zeta\omega_n s + \omega_n^2}, \quad 0 < \zeta < 1, \]

the approximate relations between the damping ratio \( \zeta \), natural frequency \( \omega_n \), and performance metrics are

\[ M_p \approx \exp\!\left( -\frac{\pi \zeta}{\sqrt{1 - \zeta^2}} \right), \quad t_s \approx \frac{4}{\zeta \omega_n}. \]

To satisfy \( M_p \leq 10\% \), it is common to select \( \zeta \approx 0.7 \), which yields a reasonably fast yet well-damped response. With \( \zeta = 0.7 \) and \( t_s \leq 1.0 \, \text{s} \), we choose

\[ \omega_n = \frac{4}{\zeta t_s} = \frac{4}{0.7 \times 1.0} \approx 5.714 \, \text{rad/s}. \]

Thus the targeted closed-loop characteristic polynomial (for a dominant second-order behavior) is

\[ s^2 + 2 \zeta \omega_n s + \omega_n^2 = s^2 + 8 s + 32.653 \quad (\text{numerically}). \]

Because \( G(s) \) already contains an integrator, the requirement \( e_{ss} = 0 \) for step input can be met using a PD controller without additional integral action.

4. PD Controller Structure and Pole-Placement Design

We select a PD controller of the form

\[ C(s) = K_p (1 + T_d s), \]

where \( K_p \) is the proportional gain and \( T_d \) is the derivative time constant. With unity feedback, the open-loop transfer function is

\[ L(s) = C(s) G(s) = \frac{K_p (1 + T_d s)}{s(0.5 s + 1)}. \]

The closed-loop transfer function is \( T(s) = \dfrac{L(s)}{1 + L(s)} \), so the characteristic equation of the closed loop is

\[ 1 + L(s) = 0 \;\Longleftrightarrow\; s(0.5 s + 1) + K_p (1 + T_d s) = 0. \]

Expanding and collecting terms:

\[ 0.5 s^2 + s + K_p T_d s + K_p = 0 \;\Longleftrightarrow\; 0.5 s^2 + (1 + K_p T_d) s + K_p = 0. \]

To match the desired second-order polynomial \( s^2 + 2 \zeta \omega_n s + \omega_n^2 \), we multiply it by \( 0.5 \) (to match the leading coefficient) and equate coefficients:

\[ 0.5 s^2 + (1 + K_p T_d) s + K_p \equiv 0.5 s^2 + (2 \zeta \omega_n \cdot 0.5) s + 0.5 \omega_n^2. \]

Therefore,

\[ 1 + K_p T_d = \zeta \omega_n, \quad K_p = 0.5 \omega_n^2. \]

Solving for \( K_p \) and \( T_d \):

\[ K_p = 0.5 \omega_n^2, \quad T_d = \frac{\zeta \omega_n - 1}{K_p}. \]

Using the numerical value \( \omega_n \approx 5.714 \) and \( \zeta = 0.7 \), we obtain

\[ K_p = 0.5 \times (5.714)^2 \approx 16.33, \quad T_d = \frac{0.7 \times 5.714 - 1}{16.33} \approx 0.184 \, \text{s}. \]

These gains place the closed-loop poles at approximately \( s = -4 \pm j4.08 \), satisfying the specified settling time and overshoot constraints under the second-order approximation.

5. Steady-State Error and Disturbance Considerations

For a unity-feedback system with open-loop transfer function \( L(s) = C(s) G(s) \), the error for a unit step reference is

\[ E(s) = \frac{1}{s} \cdot \frac{1}{1 + L(s)}, \quad e_{ss} = \lim_{t \to \infty} e(t) = \lim_{s \to 0} s E(s) = \lim_{s \to 0} \frac{1}{1 + L(s)}. \]

With \( C(s) = K_p (1 + T_d s) \) and \( G(s) = \dfrac{1}{s(0.5 s + 1)} \), we have

\[ L(s) = \frac{K_p (1 + T_d s)}{s(0.5 s + 1)} \sim \frac{K_p}{s} \quad \text{as } s \to 0, \]

so \( L(s) \to \infty \) and \( e_{ss} = 0 \) for step input. The system remains type 1 (one integrator), so it has finite nonzero error for ramp inputs, which is acceptable for a position tracking application with step-like references.

For a constant load torque modeled as an additive disturbance at the plant input, the low-frequency sensitivity function \( S(s) \) governs disturbance rejection:

\[ S(s) = \frac{1}{1 + L(s)}, \quad T(s) = \frac{L(s)}{1 + L(s)}. \]

Our design yields large loop gain around DC (due to the integrator), so \( |S(j\omega)| \) is small at low frequencies, implying good rejection of slowly-varying torque disturbances.

6. Frequency-Domain Validation

Frequency-domain validation is carried out by examining the Bode plot of \( L(s) \) and the magnitude of \( S(j\omega) \) and \( T(j\omega) \). Key quantities:

  • Gain crossover frequency \( \omega_{gc} \) such that \( |L(j\omega_{gc})| = 1 \).
  • Phase margin \( \phi_m \) at \( \omega_{gc} \).
  • Bandwidth of \( T(j\omega) \), approximately the closed-loop response speed.

Because the closed-loop poles are at approximately \( -4 \pm j4.08 \), the dominant dynamics are around \( \omega_d \approx 4.08 \). In a well-designed loop, \( \omega_{gc} \) lies on the order of \( \omega_n \), which is consistent with our pole placement.

In practice, software tools (Chapter 28) are used to compute gain and phase margins numerically and verify that the controller provides sufficient robustness to model uncertainty and unmodeled high-frequency dynamics.

7. Python Implementation (Using python-control)

We now prototype the design in Python using the control library (also known as python-control), together with numpy and matplotlib for numerical computations and plotting.


import numpy as np
import matplotlib.pyplot as plt
import control as ctrl  # python-control

# Plant parameters and specs
T = 0.5
K = 1.0
zeta = 0.7
ts_des = 1.0
wn = 4.0 / (zeta * ts_des)   # desired natural frequency

Kp = 0.5 * wn**2
Td = (zeta * wn - 1.0) / Kp

print(f"Kp = {Kp:.3f}, Td = {Td:.3f}, wn = {wn:.3f}")

# Define transfer functions
s = ctrl.TransferFunction.s
G = 1 / (s * (T * s + 1))     # normalized DC motor position plant
C = Kp * (1 + Td * s)         # PD controller

L = C * G
Tcl = ctrl.feedback(L, 1)     # closed-loop from reference to theta

# Step response
t = np.linspace(0.0, 5.0, 1000)
t_out, y_out = ctrl.step_response(Tcl, t)

plt.figure()
plt.plot(t_out, y_out)
plt.xlabel("Time [s]")
plt.ylabel("Position response theta(t)")
plt.title("Closed-loop step response (Python prototype)")
plt.grid(True)

# Sensitivity and complementary sensitivity
S = ctrl.feedback(1, L)
Tmat = ctrl.feedback(L, 1)

w = np.logspace(-2, 2, 500)
magS, phaseS, wS = ctrl.bode(S, w, Plot=False)
magT, phaseT, wT = ctrl.bode(Tmat, w, Plot=False)

plt.figure()
plt.semilogx(wS, 20 * np.log10(magS), label="S(jw)")
plt.semilogx(wT, 20 * np.log10(magT), label="T(jw)")
plt.xlabel("Frequency [rad/s]")
plt.ylabel("Magnitude [dB]")
plt.title("Sensitivity and complementary sensitivity")
plt.grid(True, which="both")
plt.legend()

plt.show()
      

This Python prototype confirms that the PD gains derived analytically deliver the desired step response and that the sensitivity function has a low magnitude near DC, indicating good disturbance rejection.

flowchart TD
  M["Model G(s), C(s)"] --> PY["Python prototype (python-control)"]
  PY --> VAL["Check specs: ts, Mp, ess, margins"]
  VAL --> CODE["Export gains to C++ / Java / MATLAB / Mathematica"]
        

8. C++ Implementation (Numerical Simulation Loop)

For embedded or real-time applications, C++ is frequently used with a numerical integration scheme to simulate or implement the controller. Using the time-domain model corresponding to \( G(s) = \dfrac{1}{s(0.5 s + 1)} \), we obtain

\[ 0.5 \ddot{\theta}(t) + \dot{\theta}(t) = u(t). \]

Setting \( x_1 = \theta \), \( x_2 = \dot{\theta} \), and using the PD controller with a constant step reference \( r(t) = 1 \):

\[ u(t) = K_p \left(r(t) - x_1(t) - T_d x_2(t)\right), \]

the state-space form is

\[ \dot{x}_1 = x_2, \quad \dot{x}_2 = \frac{-x_2 + u}{0.5}. \]

A simple explicit Euler discretization with step size \( h \) yields:


#include <iostream>
#include <vector>

int main() {
    const double T = 0.5;
    const double zeta = 0.7;
    const double ts_des = 1.0;
    const double wn = 4.0 / (zeta * ts_des);

    const double Kp = 0.5 * wn * wn;
    const double Td = (zeta * wn - 1.0) / Kp;

    std::cout << "Kp = " << Kp << ", Td = " << Td << std::endl;

    double h = 1e-3;     // integration time step
    double t_end = 5.0;
    int N = static_cast<int>(t_end / h);

    // state: x1 = theta, x2 = theta_dot
    double x1 = 0.0;
    double x2 = 0.0;
    double r = 1.0;      // unit step reference

    std::vector<double> t_vec;
    std::vector<double> y_vec;

    t_vec.reserve(N + 1);
    y_vec.reserve(N + 1);

    t_vec.push_back(0.0);
    y_vec.push_back(x1);

    for (int k = 0; k < N; ++k) {
        double t = (k + 1) * h;

        // PD control law: u = Kp * (e - Td * theta_dot)
        double e = r - x1;
        double u = Kp * (e - Td * x2);

        // continuous-time dynamics
        double x1_dot = x2;
        double x2_dot = (-x2 + u) / T;

        // Euler integration
        x1 += h * x1_dot;
        x2 += h * x2_dot;

        t_vec.push_back(t);
        y_vec.push_back(x1);
    }

    // Print a few samples or export to file
    for (size_t k = 0; k < t_vec.size(); k += 500) {
        std::cout << t_vec[k] << "  " << y_vec[k] << std::endl;
    }

    return 0;
}
      

Here, C++ serves as a stand-in for an embedded implementation: the same control law can be executed on a microcontroller, with the Euler step replaced by sampling and hold of the real plant.

9. Java Implementation (Simulation and Logging)

Java is often used in higher-level supervisory or simulation tools, where we may still need a basic numerical integration of the same ODEs. The controller law and plant dynamics mirror the C++ version.


public class PdMotorSimulation {
    public static void main(String[] args) {
        double T = 0.5;
        double zeta = 0.7;
        double tsDes = 1.0;
        double wn = 4.0 / (zeta * tsDes);

        double Kp = 0.5 * wn * wn;
        double Td = (zeta * wn - 1.0) / Kp;

        System.out.printf("Kp = %.3f, Td = %.3f, wn = %.3f%n", Kp, Td, wn);

        double h = 1e-3;
        double tEnd = 5.0;
        int N = (int) (tEnd / h);

        double x1 = 0.0;  // theta
        double x2 = 0.0;  // theta_dot
        double r = 1.0;

        for (int k = 0; k <= N; ++k) {
            double t = k * h;

            double e = r - x1;
            double u = Kp * (e - Td * x2);

            double x1Dot = x2;
            double x2Dot = (-x2 + u) / T;

            x1 += h * x1Dot;
            x2 += h * x2Dot;

            if (k % 1000 == 0) {
                System.out.printf("t = %.3f, theta = %.4f%n", t, x1);
            }
        }
    }
}
      

This Java simulation can be integrated into a GUI-based tool or used for regression tests that verify closed-loop behavior against reference data generated in Python or MATLAB.

10. MATLAB/Simulink Implementation

MATLAB is a standard environment for linear control design. We directly encode the transfer functions and PD controller, and we can represent the same loop in Simulink using standard blocks (Transfer Fcn, Sum, Gain).


% Plant and design parameters
T = 0.5;
zeta = 0.7;
ts_des = 1.0;
wn = 4 / (zeta * ts_des);

Kp = 0.5 * wn^2;
Td = (zeta * wn - 1) / Kp;

s = tf('s');

G = 1 / (s * (T*s + 1));     % DC motor position plant
C = Kp * (1 + Td * s);       % PD controller

L = C * G;
Tcl = feedback(L, 1);

% Step response
figure;
step(Tcl);
grid on;
title('Closed-loop step response (MATLAB)');

% Sensitivity and complementary sensitivity
S = feedback(1, L);
Tmat = feedback(L, 1);

figure;
bodemag(S, Tmat);
grid on;
legend('S(s)', 'T(s)');
title('Sensitivity and complementary sensitivity');

% Simulink hint:
% - Create a new model.
% - Add a Step block for reference r(t).
% - Add Sum block for error e(t) = r - theta.
% - Implement C(s) as a Transfer Fcn or as series of Gain and Derivative blocks.
% - Add Transfer Fcn block for G(s).
% - Close the loop with a second Sum block.
% - Use Scope block to visualize theta(t) and e(t).
      

In Simulink, the continuous-time implementation of the PD controller can use a filtered derivative for robustness, replacing \( T_d s \) by \( \dfrac{T_d s}{\alpha T_d s + 1} \) with \( \alpha \in (0,1) \), as covered in the PID implementation lessons.

11. Wolfram Mathematica Implementation

Wolfram Mathematica offers symbolic and numeric tools for linear control, for example via TransferFunctionModel and FeedbackConnect. Below we reproduce the PD design and plot the step response.


(* Parameters *)
T = 0.5;
zeta = 0.7;
tsDes = 1.0;
wn = 4.0/(zeta*tsDes);

Kp = 0.5*wn^2;
Td = (zeta*wn - 1.0)/Kp;

(* Laplace variable s *)
s = LaplaceTransformVariable;

(* Plant and controller *)
G = TransferFunctionModel[1/(s*(T*s + 1)), s];
C = TransferFunctionModel[Kp*(1 + Td*s), s];

(* Closed-loop transfer function from reference to theta *)
Tsys = FeedbackConnect[C*G, 1];

(* Step response *)
StepResponsePlot[Tsys, {t, 0, 5},
  PlotLabel -> "Closed-loop step response (Mathematica)",
  Frame -> True,
  GridLines -> Automatic
]

(* Poles of closed loop *)
Poles[Tsys]
      

Mathematica is particularly useful for symbolic manipulation of characteristic polynomials, allowing formal verification of the relations between controller gains and desired closed-loop poles derived earlier.

12. Problems and Solutions

Problem 1 (PD Design Formula Derivation): For the plant \( G(s) = \dfrac{1}{s(0.5 s + 1)} \) and PD controller \( C(s) = K_p (1 + T_d s) \), derive the formulas relating \( K_p \) and \( T_d \) to the desired second-order parameters \( \zeta \) and \( \omega_n \).

Solution: The closed-loop characteristic equation is

\[ 1 + L(s) = 0 \;\Longleftrightarrow\; s(0.5 s + 1) + K_p (1 + T_d s) = 0, \]

which expands to

\[ 0.5 s^2 + s + K_p T_d s + K_p = 0 \;\Longleftrightarrow\; 0.5 s^2 + (1 + K_p T_d) s + K_p = 0. \]

We match this with \( 0.5(s^2 + 2\zeta\omega_n s + \omega_n^2) \), giving

\[ 1 + K_p T_d = \zeta \omega_n, \quad K_p = 0.5 \omega_n^2. \]

Solving for \( T_d \):

\[ T_d = \frac{\zeta \omega_n - 1}{K_p}. \]

Problem 2 (Steady-State Error Classification): Show that the PD-controlled system remains type 1 and therefore has zero steady-state error to step input, and nonzero finite steady-state error to a ramp input.

Solution: The open-loop transfer function is \( L(s) = C(s) G(s) = \dfrac{K_p (1 + T_d s)}{s(0.5 s + 1)} \). As \( s \to 0 \), we have \( L(s) \sim \dfrac{K_p}{s} \), hence \( \lim_{s \to 0} L(s) = \infty \). The system has one pole at the origin and no zeros cancelling it, so it is type 1. For a unit step,

\[ e_{ss} = \lim_{s \to 0} \frac{1}{1 + L(s)} = 0. \]

For a unit ramp input \( R(s) = 1/s^2 \), the error is

\[ E(s) = \frac{1}{s^2} \cdot \frac{1}{1 + L(s)}, \quad e_{ss} = \lim_{s \to 0} s E(s) = \lim_{s \to 0} \frac{1}{s(1 + L(s))}. \]

Since \( L(s) \sim K_p / s \), we obtain a finite nonzero limit \( e_{ss} = 1/K_v \), where \( K_v = \lim_{s \to 0} s L(s) = K_p \). Thus the ramp error is finite and nonzero, consistent with type-1 behavior.

Problem 3 (Constraint on \( T_d \) for Positivity): Show that the derivative time \( T_d \) is positive if and only if \( \zeta \omega_n > 1 \). Interpret this condition in terms of achievable damping for the given plant.

Solution: From the design formula \( T_d = (\zeta \omega_n - 1)/K_p \), and since \( K_p > 0 \), we require \( \zeta \omega_n - 1 > 0 \), i.e.

\[ \zeta \omega_n > 1. \]

This inequality states that for a given plant (leading coefficient fixed at \( 0.5 \)), only sufficiently fast and/or well-damped closed-loop pole locations admit a positive derivative term. For slower or less damped choices, a pure proportional controller already saturates the achievable damping and a positive physical derivative action cannot be used to match the desired polynomial exactly.

Problem 4 (Effect of Plant Time Constant Variation): Suppose that the true plant is \( \tilde{G}(s) = \dfrac{1}{s(\tilde{T} s + 1)} \), with \( \tilde{T} = 0.6 \) instead of \( T = 0.5 \), but the controller gains are still designed using \( T = 0.5 \). Qualitatively, does the closed-loop become slower or faster, and is the damping ratio increased or decreased?

Solution: The characteristic polynomial with uncertainty becomes

\[ \tilde{T} s^2 + (1 + K_p T_d) s + K_p = 0, \]

while the gains \( K_p, T_d \) were computed to match the coefficients for \( T = 0.5 \). Increasing \( \tilde{T} \) increases the inertia term in the denominator, reducing the effective natural frequency and hence slowing down the response. The ratio between the linear and quadratic coefficients changes in such a way that the damping ratio is typically reduced (the system tends to become more oscillatory), illustrating sensitivity of pole placement to modeling errors in \( T \).

Problem 5 (Controller Implementation Order): Consider a digital implementation of the PD controller with sampling period \( h = 2 \, \text{ms} \). Discuss the consequences of discretizing the plant and controller versus using a continuous-time model for design and then applying a straightforward backward difference approximation for the derivative term. Why is it important to check the discrete-time closed-loop poles?

Solution: Designing in continuous time and then discretizing may shift the effective pole locations because the mapping from the \( s \)-plane to the \( z \)-plane is nonlinear. Backward difference approximations introduce phase lag and gain variations at higher frequencies; if these are not accounted for, the actual discrete system can have reduced phase margin and potentially unstable behavior. Therefore, after discretization, one should form the discrete-time transfer functions \( G(z) \) and \( C(z) \), compute the closed-loop poles in the \( z \)-plane, and verify that they lie inside the unit disk with sufficient damping, consistent with the original continuous-time specifications.

13. Summary

In this mini-project, we completed an end-to-end linear controller design for a DC-motor-driven rotary stage. Starting from mechanical and electrical equations, we obtained a normalized plant model, translated time-domain specifications into desired pole locations, and used polynomial matching to compute PD gains analytically. We then verified steady-state error and disturbance rejection properties using sensitivity functions, and implemented the design in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica.

The workflow exemplifies how classical linear control theory links modeling, analysis, and implementation, and how analytical formulas derived in earlier chapters are used in a realistic design task.

14. References

  1. Nyquist, H. (1932). Regeneration theory. Bell System Technical Journal, 11(1), 126–147.
  2. Bode, H.W. (1945). Network Analysis and Feedback Amplifier Design. Van Nostrand, New York.
  3. Evans, W.R. (1948). Graphical analysis of control systems. Transactions of the AIEE, 67(1), 547–551.
  4. Ziegler, J.G., & Nichols, N.B. (1942). Optimum settings for automatic controllers. Transactions of the ASME, 64(11), 759–768.
  5. Freudenberg, J.S., & Looze, D.P. (1985). Right half plane poles and zeros and design tradeoffs in feedback systems. IEEE Transactions on Automatic Control, 30(6), 555–565.
  6. Middleton, R.H., & Goodwin, G.C. (1990). Digital control and estimation: A unified approach. Prentice Hall.
  7. Doyle, J.C., Francis, B.A., & Tannenbaum, A.R. (1992). Feedback Control Theory. Macmillan.
  8. Åström, K.J., & Hägglund, T. (1995). PID Controllers: Theory, Design, and Tuning. Instrument Society of America.
  9. Skogestad, S., & Postlethwaite, I. (1996). Multivariable Feedback Control: Analysis and Design. Wiley.
  10. Kuo, B.C., & Golnaraghi, F. (2002). Automatic Control Systems. Wiley.