SISO point control (regulation)

This notebook illustrates examples of single-input, single-output point control. All user input sections are highlighted as cell boxes.

Definitions

Generic definitions

Turn off annoying warnings messages

Off[General::spell1];
Off[General::spell];

Load in necessary packages

<< Graphics`Colors`;

Differential equation to characteristic equation

DiffEqToCharEq[
          eq_] := 
eq /.  {Derivative[i_][x][
                t] -> s^i, x[t] -> 1};

Display DSolve output

DisplayDSolve[eqs_] := 
      x == (x  /.  DSolve[eqs, x, t])[[1, 2]]  //  
        TraditionalForm

Process DSolve output

ProcessDSolve[
        eqs_] := 
(x  /.  
              Flatten[DSolve[eqs, x, t]])[[2]]  //  
        ExpandAll

Root locus plotting

One root locus for one root

OneRootLocusOneRoot[root_, const_, vary_, 
          opts___] := 
Module[{tmp, 
            fnc}, 
tmp = 
            root  /.  const[[1]] -> 
                const[[2]]; 
fnc  =  {Re[tmp], 
              Im[tmp]}; 
Show[
            ParametricPlot[Evaluate[fnc], vary, opts], 
            Show[Graphics[{{Blue, PointSize[ .02], 
                    Point[{Re[tmp], Im[tmp]} /. 
                        vary[[1]] -> vary[[2]]]}, {Magenta, 
                    PointSize[ .02], 
                    Point[{Re[tmp], Im[tmp]} /. 
                        vary[[1]] -> vary[[3]]]}}], 
              DisplayFunction -> 
                Identity]]
];

Corresponding gain line

OneGainLineOneRoot[root_, const_, vary_, 
          opts___] := 
 {{opts, 
              Line[{{const[[2]], vary[[2]]}, {const[[2]], 
                    vary[[3]]}}]}, {Blue, PointSize[ .02], 
              Point[{const[[2]], vary[[2]]}]}, {Magenta, 
              PointSize[ .02], 
              Point[{const[[2]], vary[[3]]}]}}  //  N;

All root loci for one root

RootLociOneRoot[root_, v1_, v2_, opts___] := 
Table[
          OneRootLocusOneRoot[
            root, {v1[[1]], i}, {v2[[1]], v2[[2]], 
              v2[[3]]},  Frame -> True,  
            DisplayFunction -> Identity, opts], {i, v1[[2]], 
            v1[[3]], v1[[4]]}];

Corresponding gain lines

GainLinesOneRoot[root_, v1_, v2_, 
          opts___] := 
Table[
          OneGainLineOneRoot[
            root, {v1[[1]], i}, {v2[[1]], v2[[2]], v2[[3]]}, 
            opts], {i, v1[[2]], v1[[3]], v1[[4]]}];

Root loci for both roots

RootLoci[roots_, v1_, v2_,  opts___] := 
Show[
          RootLociOneRoot[roots[[1]], v1, v2, PlotStyle -> Red]~Join~
            RootLociOneRoot[roots[[2]], v1, v2, PlotStyle -> Green], 
          Axes -> False, PlotLabel ->  "Complex Plane", 
          opts];

Corresponding gain lines

GainLines[roots_, v1_, v2_,  opts___] := 
Show[
          Graphics[GainLinesOneRoot[roots[[1]], v1, v2, Black]],  
          AspectRatio -> Automatic,  Frame -> True, 
          PlotLabel -> "Gain trajectories", opts];

Uncontrolled second-order linear system

Differential equation

Define simple second-order differential equation

diffeq  =  
        m  x''[t]  +  b  x'[t] +  k  x[t]  ==  
          0;

Display differential equation

diffeq  //  TraditionalForm

Differential equation solution

Solution to differential equation

DisplayDSolve[diffeq]

Set up initial conditions

initCond  =  {x[0] == x0, x'[0] == 
            dx0};
\(\(displayRule\  = \ {x0 -> x\_0, 
          dx0 -> x\& . \_0};\)\)

Solution to differential equation with initial conditions

DisplayDSolve[{diffeq, initCond}]  /.  displayRule

Characteristic equation

Compute characteristic equation

chareq  =  DiffEqToCharEq[diffeq];

Display characteristic equation

chareq  //  TraditionalForm

Roots of characteristic equation

roots  =  (s /. #) &  /@  Solve[chareq, s];

Display roots of characteristic equation

roots  //  TraditionalForm

Solution classes

Characteristic equation normalized with respect to m

chareq  =  DiffEqToCharEq[diffeq];
chareq1  = (chareq[[1]]/m  ==  0  //  
            ExpandAll);
chareq1  //  TraditionalForm

Equivalent representation (ζ = damping ration, ω = natural frequency)

chareq2  =  
        s^2  +  2  ζ ω s +  ω^2  ==  
          0;
chareq2  //  TraditionalForm

Roots of above characteristic equation

(s  /.  #) &  /@  Solve[chareq2, s]  //  
      TraditionalForm

Solution classes

\(\(\({{ζ < 1, \ "underdamped"}, {ζ == 
            1, "critically damped"}, {ζ > 
            1, "overdamped"}}\)\(\ \)\(//\)\(\ \)\(TableForm\)\(\ \)\)\)

Solution for ζ and ω in terms of system parameters (m,  b,  k)

ωs  = Solve[(ω^2) ==  k/m, ω]  //  
            Flatten  //  Last;
ωs  //  TraditionalForm

ζs = Solve[
              2  ζ ω  ==  
                  b/m  /.  ωs, ζ]  //  Flatten  //  
          Last;
ζs  // TraditionalForm

Overdampled example (unequal negative real roots, ζ > 1)

Numeric values

vals  =  {m -> 1,  b -> 9,  k -> 6,  x0 -> 1,  
            dx0 -> 0}  //  N;

Damping ratio and natural frequency

{ζ, ω}  /.  {ζs, ωs}  /.  
      vals

Corresponding roots

roots  /.  vals

Differential equation solution

fx = 
        ProcessDSolve[{diffeq, initCond}  /.  
            vals];
x[t]  == fx  //  TraditionalForm

Plot length

tmax  =  10. ;

Plot solution

g1 = 
        Plot[fx, {t, 0, tmax},  PlotStyle -> Blue, Frame -> True, 
          PlotRange -> All];

Underdamped system (conjugate complex roots with negative real parts, ζ < 1)

Numeric values

vals  =  {m -> 1,  b -> 1,  k -> 4,  x0 -> 1,  
            dx0 -> 0}  //  N;

Damping ratio and natural frequency

{ζ, ω}  /.  {ζs, ωs}  /.  
      vals

Corresponding roots

roots  /.  vals

Differential equation solution

fx = 
        ProcessDSolve[{diffeq, initCond}  /.  
            vals];
x[t]  == fx  //  TraditionalForm

Plot length

tmax  =  10. ;

Plot solution

g2 = 
        Plot[fx, {t, 0, tmax},  PlotStyle -> Red, Frame -> True, 
          PlotRange -> All];

Critically damped system (equal negative real roots, ζ = 1)

Numeric values

The values below are exact, because of numeric instability in DSolve[ ]. Use NDSolve[ ] for non-exact values.

vals  =  {m -> 1,  b -> 4,  k -> 4,  x0 -> 1,  
          dx0 -> 0};

Damping ratio and natural frequency

{ζ, ω}  /.  {ζs, ωs}  /.  
      vals

Corresponding roots

roots  /.  vals

Differential equation solution

fx = 
        ProcessDSolve[{diffeq, initCond}  /.  
            vals];
x[t]  == fx  //  TraditionalForm

Plot length

tmax  =  10. ;

Plot solution

g3 = 
        Plot[fx, {t, 0, tmax},  PlotStyle -> Magenta, 
          Frame -> True, PlotRange -> All];

Compare three stable examples

Show[g1, g2, g3];

Unstable system (roots have positive real parts)

Numeric values

The values below are exact, because of numeric instability in DSolve[ ]. Use NDSolve[ ] for non-exact values.

vals  =  {m -> 1,  b -> -1,  k -> 10,  
          x0 -> 1,  dx0 -> 0};

Damping ratio and natural frequency

{ζ, ω}  /.  {ζs, ωs}  /.  
      vals

Corresponding roots

roots  /.  vals

Differential equation solution

fx = 
        ProcessDSolve[{diffeq, initCond}  /.  
            vals];
x[t]  == fx  //  TraditionalForm

Plot length

tmax  =  10. ;

Plot solution

g4 = 
        Plot[fx, {t, 0, tmax},  PlotStyle -> Green, 
          Frame -> True];

PD control  (stable, linear system)

Differential equation

Define simple second-order differential equation with input f

diffeq  =  
        m  x''[t]  +  b  x'[t] +  k  x[t]  ==  
          f;

Display differential equation

diffeq  //  TraditionalForm

Numeric system values

The values below are exact, because of numeric instability in DSolve[ ]. Use NDSolve[ ] for non-exact values.

vals  =  {m -> 1,  b -> 1,  k -> 1,  x0 -> 1,  
          dx0 -> 0};

Zero control (f = 0)

Control input

control = {f -> 0};

Damping ratio and natural frequency

{ζ, ω}  /.  {ζs, ωs}  /.  
      vals

Corresponding roots

roots  /.  vals

Differential equation solution

fx = 
        ProcessDSolve[{diffeq, initCond}  /.  vals  /.  
              control  //  N];
x[t]  == fx  //  TraditionalForm

Plot length

tmax  =  10. ;

Plot solution

g1 = 
        Plot[fx, {t, 0, tmax},  PlotStyle -> Blue, Frame -> True, 
          PlotRange -> All];

PD control (\(TraditionalForm\`f\  = \ \(-k\_p\) x\  - \ \(k\_v\) \(x\& . \)\))

Control input

control = {f -> -kp x[t]  -  kv x'[t]};

Corresponding controlled system

controlledSystem  =  diffeq  /.  control;

Display controlled system equation

controlledSystem  //  TraditionalForm

controlledSystem  =  
        controlledSystem  //  Simplify;
controlledSystem  //  TraditionalForm

Desired controlled system characteristics

The values below are exact, because of numeric instability in DSolve[ ]. Use NDSolve[ ] for non-exact values.

desiredCharacteristics  =  {ω -> 4, ζ -> 
            1};

Computation of gain values for characteristic equation (s^2) + 2  ζsω + (ω^2) = 0

eqs  =  {ω == (ω /. ωs), ζ  ==  (ζ /. ζs)};
eqs2  =  
        eqs  /.  {k -> (k + kp), 
            b -> (b + kv)};
eqs2  //  TraditionalForm

gainsSymbolic  =  Solve[eqs2, {kp, kv}]  //  Flatten

gains  =  gainsSymbolic  /.  vals  /. 
        desiredCharacteristics

Controlled system roots

controlledCharEq  =  
        DiffEqToCharEq[controlledSystem];
controlledRoots  =  (s /. #) &  /@  
          Solve[controlledCharEq, s];

Display controlled system roots

controlledRoots  /.  
        vals /. gains

Differential equation solution

fx = 
        ProcessDSolve[{controlledSystem, initCond}  /.  vals  /.  
            gains];
x[t]  == fx  //  TraditionalForm

Plot length

tmax  =  10. ;

Plot solution

g2 = 
        Plot[fx, {t, 0, tmax},  PlotStyle -> Red, Frame -> True, 
          PlotRange -> All];

Compare to unforced system

Show[g1, g2];

Control f

dfx  =  D[fx, t];
fn  =  f  /.  control  /.  
            gains /.  {x[t] -> fx,  x'[t] -> 
              dfx};
fn  //  TraditionalForm

Plot[fn, {t, 0, tmax},  PlotStyle -> Magenta, 
        Frame -> True, PlotRange -> All];

Root locus plot

kpRange = {kp, 15, 15, 1};
kvRange = {kv, 0, 20};

g3  = 
        GainLines[controlledRoots  /.  vals, kpRange, kvRange,  
          DisplayFunction -> Identity,  
          PlotRange -> {Drop[kvRange, 1] + {-1, 1}, 
              Drop[kvRange, 1] + {-1, 1}}];
g4  =  
        RootLoci[controlledRoots  /.  vals, kpRange, kvRange,  
          DisplayFunction -> Identity,  
          PlotRange -> All];
Show[GraphicsArray[{{g3}, {g4}}], 
        DisplayFunction -> $DisplayFunction];

PD control  (unstable, linear system)

Differential equation

Define simple second-order differential equation with input f

diffeq  =  
        m  x''[t]  +  b  x'[t] +  k  x[t]  ==  
          f;

Display differential equation

diffeq  //  TraditionalForm

Numeric system values

The values below are exact, because of numeric instability in DSolve[ ]. Use NDSolve[ ] for non-exact values.

vals  =  {m -> 1,  b -> -1,  k -> 10,  
          x0 -> 1,  dx0 -> 0};

Zero control (f = 0)

Control input

control = {f -> 0};

Damping ratio and natural frequency

{ζ, ω}  /.  {ζs, ωs}  /.  
      vals

Corresponding roots

roots  /.  vals

Differential equation solution

fx = 
        ProcessDSolve[{diffeq, initCond}  /.  vals  /.  
              control  //  N];
x[t]  == fx  //  TraditionalForm

Plot length

tmax  =  10. ;

Plot solution

g1 = 
        Plot[fx, {t, 0, tmax},  PlotStyle -> Blue, Frame -> True, 
          PlotRange -> All];

PD control (\(TraditionalForm\`f\  = \ \(-k\_p\) x\  - \ \(k\_v\) \(x\& . \)\))

Control input

control = {f -> -kp x[t]  -  kv x'[t]};

Corresponding controlled system

controlledSystem  =  diffeq  /.  control;

Display controlled system equation

controlledSystem  //  TraditionalForm

controlledSystem  =  
        controlledSystem  //  Simplify;
controlledSystem  //  TraditionalForm

Desired controlled system characteristics

The values below are exact, because of numeric instability in DSolve[ ]. Use NDSolve[ ] for non-exact values.

desiredCharacteristics  =  {ω -> 4, ζ -> 
            1};

Computation of gain values for characteristic equation (s^2) + 2  ζsω + (ω^2) = 0

eqs  =  {ω == (ω /. ωs), ζ  ==  (ζ /. ζs)};
eqs2  =  
        eqs  /.   {k -> 
              Coefficient[Expand[controlledSystem[[1]]], x[t]], 
            b -> 
              Coefficient[
                Expand[controlledSystem[[1]]], x'[
                  t]]};
eqs2  //  PowerExpand  //  TraditionalForm

gainsSymbolic  =  Solve[eqs2, {kp, kv}]  //  Flatten

gains  =  gainsSymbolic  /.  vals  /. 
        desiredCharacteristics

Controlled system roots

controlledCharEq  =  
        DiffEqToCharEq[controlledSystem];
controlledRoots  =  (s /. #) &  /@  
          Solve[controlledCharEq, s];

Display controlled system roots

controlledRoots  /.  
        vals /. gains

Differential equation solution

fx = 
        ProcessDSolve[{controlledSystem, initCond}  /.  vals  /.  
            gains];
x[t]  == fx  //  TraditionalForm

Plot length

tmax  =  10. ;

Plot solution

g2 = 
        Plot[fx, {t, 0, tmax},  PlotStyle -> Red, Frame -> True, 
          PlotRange -> All];

Compare to unforced system

Show[g1];

Control f

dfx  =  D[fx, t];
fn  =  f  /.  control  /.  
            gains /.  {x[t] -> fx,  x'[t] -> 
              dfx};
fn  //  TraditionalForm

Plot[fn, {t, 0, tmax},  PlotStyle -> Magenta, 
        Frame -> True, PlotRange -> All];

Root locus plot

kpRange = {kp, 6, 6, 1};
kvRange = {kv, 0, 10};

g3  = 
        GainLines[controlledRoots  /.  vals, kpRange, kvRange,  
          DisplayFunction -> Identity,  
          PlotRange -> {Drop[kvRange, 1] + {-1, 1}, 
              Drop[kvRange, 1] + {-1, 1}}];
g4  =  
        RootLoci[controlledRoots  /.  vals, kpRange, kvRange,  
          DisplayFunction -> Identity,  
          PlotRange -> All];
Show[GraphicsArray[{{g3}, {g4}}], 
        DisplayFunction -> $DisplayFunction];

Partitioned control (stable, linear system)

Differential equation

Define simple second-order differential equation with input f

diffeq  =  
        m  x''[t]  +  b  x'[t] +  k  x[t]  ==  
          f;

Display differential equation

diffeq  //  TraditionalForm

Numeric system values

The values below are exact, because of numeric instability in DSolve[ ]. Use NDSolve[ ] for non-exact values.

vals  =  {m -> 1,  b -> 1,  k -> 1,  x0 -> 1,  
          dx0 -> 0};

Zero control (f = 0)

Control input

control = {f -> 0};

Damping ratio and natural frequency

{ζ, ω}  /.  {ζs, ωs}  /.  
      vals

Corresponding roots

roots  /.  vals

Differential equation solution

fx = 
        ProcessDSolve[{diffeq, initCond}  /.  vals  /.  
              control  //  N];
x[t]  == fx  //  TraditionalForm

Plot length

tmax  =  10. ;

Plot solution

g1 = 
        Plot[fx, {t, 0, tmax},  PlotStyle -> Blue, Frame -> True, 
          PlotRange -> All];

Partitioned control

Partitioned control law

Partitioned control law: 
	\(TraditionalForm\`f = \(\(α\ \(f\&^\)\)\(\ \)\(+\)\(β\)\(\ \)\)\) where,
	
	α = m 
	\(TraditionalForm\`β\  = \ b\ \(x\& . \)\  + \ k\ x\) 
	\(TraditionalForm\`\(f\&^\)\  = \ \(-k\_p\) 
            x\  - \ \(k\_v\) \(x\& . \)\)

control = {f -> α fp  +  β , 
          fp -> -kp x[t]  -  kv x'[t],  α -> 
            m, β -> b x'[t]  +  k x[t]};

Corresponding controlled system

controlledSystem  =  diffeq  //.  control;

Display controlled system equation

controlledSystem  //  TraditionalForm

controlledSystem  =  (controlledSystem  //  
            Simplify);
controlledSystem  //  TraditionalForm

Desired controlled system characteristics

The values below are exact, because of numeric instability in DSolve[ ]. Use NDSolve[ ] for non-exact values.

desiredCharacteristics  =  {ω -> 4, ζ -> 
            1};

Computation of gain values for characteristic equation (s^2) + 2  ζsω + (ω^2) = 0

eqs  =  {ω == (ω /. ωs), ζ  ==  (ζ /. ζs)};
eqs2  =  
        eqs  /.   {k -> 
              Coefficient[Expand[controlledSystem[[1]]], x[t]], 
            b -> 
              Coefficient[
                Expand[controlledSystem[[1]]], x'[
                  t]]};
eqs2  //  PowerExpand  //  TraditionalForm

gainsSymbolic  =  Solve[eqs2, {kp, kv}]  //  Flatten

gains  =  gainsSymbolic  /.  vals  /. 
        desiredCharacteristics

Controlled system roots

controlledCharEq  =  
        DiffEqToCharEq[controlledSystem];
controlledRoots  =  (s /. #) &  /@  
          Solve[controlledCharEq, s];

Display controlled system roots

controlledRoots  /.  
        vals /. gains

Differential equation solution

fx = 
        ProcessDSolve[{controlledSystem, initCond}  /.  vals  /.  
            gains];
x[t]  == fx  //  TraditionalForm

Plot length

tmax  =  10. ;

Plot solution

g2 = 
        Plot[fx, {t, 0, tmax},  PlotStyle -> Red, Frame -> True, 
          PlotRange -> All];

Compare to unforced system

Show[g1, g2];

Control f

dfx  =  D[fx, t];
fn  =  f  //.  control  /. vals  /.  
            gains /.  {x[t] -> fx,  x'[t] -> 
              dfx};
fn  //  TraditionalForm

Plot[fn, {t, 0, tmax},  PlotStyle -> Magenta, 
        Frame -> True, PlotRange -> All];

Partitioned control (stable, nonlinear system)

Differential equation

Define simple second-order differential equation with input f

diffeq  =  
        m  x''[t]  +  b  x'[t] +  q  x[t]^3  ==  
          f;

Display differential equation

diffeq  //  TraditionalForm

Numeric system values

vals  =  {m -> 1,  b -> 1,  q -> 1,  x0 -> 1,  
          dx0 -> 0};

Zero control (f = 0)

Control input

control = {f -> 0};

Simulation length

tmax  =  10. ;

Numerical differential equation solution

fx = 
      x  /.  (NDSolve[{diffeq, initCond}  /.  vals  /.  
                  control  //  N, x, {t, 0, tmax}]  //  
            Flatten)

Plot solution

g1 = 
        Plot[fx[t], {t, 0, tmax},  PlotStyle -> Blue, 
          Frame -> True, PlotRange -> {All, {0, 1}}];

Partitioned control

Partitioned control law

Partitioned control law: 
	\(TraditionalForm\`f = \(\(α\ \(f\&^\)\)\(\ \)\(+\)\(β\)\(\ \)\)\) where,
	
	α = m 
	\(TraditionalForm\`β\  = \ b\ \(x\& . \)\  + \ q\ x\^3\) 
	\(TraditionalForm\`\(f\&^\)\  = \ \(-k\_p\) 
            x\  - \ \(k\_v\) \(x\& . \)\)

control = {f -> α fp  +  β , 
          fp -> -kp x[t]  -  kv x'[t],  α -> 
            m, β -> b x'[t]  +  q x[t]^3};

Corresponding controlled system

controlledSystem  =  diffeq  //.  control;

Display controlled system equation

controlledSystem  //  TraditionalForm

controlledSystem  =  (controlledSystem  //  
            Simplify);
controlledSystem  //  TraditionalForm

Desired controlled system characteristics

desiredCharacteristics  =  {ω -> 4, ζ -> 
            1};

Computation of gain values for characteristic equation (s^2) + 2  ζsω + (ω^2) = 0

eqs  =  {ω == (ω /. ωs), ζ  ==  (ζ /. ζs)};
eqs2  =  
        eqs  /.   {k -> 
              Coefficient[Expand[controlledSystem[[1]]], x[t]], 
            b -> 
              Coefficient[
                Expand[controlledSystem[[1]]], x'[
                  t]]};
eqs2  //  PowerExpand  //  TraditionalForm

gainsSymbolic  =  Solve[eqs2, {kp, kv}]  //  Flatten

gains  =  gainsSymbolic  /.  vals  /. 
        desiredCharacteristics

Numerical differential equation solution

fx = 
        x  /.  (NDSolve[{controlledSystem, initCond}  /.  
                      vals  /.  gains  //  N, x, {t, 0, tmax}]  //  
              Flatten);

Plot solution

g2 = 
        Plot[fx[t], {t, 0, tmax},  PlotStyle -> Red, 
          Frame -> True, PlotRange -> All];

Compare to unforced system

Show[g2, g1];

Control f

dfx  =  D[fx[t], t];
fn  =  f  //.  control  /. vals  /.  
            gains /.  {x[t] -> fx[t],  x'[t] -> 
              dfx};

Plot[Evaluate[fn], {t, 0, tmax},  PlotStyle -> Magenta, 
        Frame -> True, PlotRange -> All];