To derive the model of the system, first, we apply Kirchhoff's voltage law around the circuit. The back EMF is proportional to the angular velocity of the motor shaft. Using the $\mathrm{DiffEquation}$ command we get:

$$\begin{gathered} \mathrm{sysElectrical} ≔ \mathrm{DiffEquation}\left(\left[L\cdot i\'\left(t\right)\+R\cdot i\left(t\right) \= V\left(t\right)-K\cdot \mathrm{\θ}\'\left(t\right)\right]\, \left[V\left(t\right)\right]\, \left[\mathrm{\θ}\left(t\right)\right]\right)\: \\ \mathrm{PrintSystem}\left(\mathrm{sysElectrical}\right) \end{gathered}$$

Second, we model the mechanical component of the motor as a second order differential system with moment of inertia $J$ and damping ratio of $b$. The torque is proportional to the current $i\(t\)$ driving the motor.

$$\begin{gathered} \mathrm{sysMechanical} ≔ \mathrm{DiffEquation}\left(\left[J\cdot \mathrm{\θ}\'\'\left(t\right)\+b\cdot \mathrm{\θ}\'\left(t\right) \= K\cdot i\left(t\right)\right]\, \left[i\left(t\right)\right]\, \left[\mathrm{\θ}\left(t\right)\right]\right)\: \\ \mathrm{PrintSystem}\left(\mathrm{sysMechanical}\right) \end{gathered}$$

The overall differential equation model in terms of the voltage, $V\left(t\right)$, and angular position, $\mathrm{\Θ}\left(t\right)$, is:

$$\begin{gathered} \mathrm{sysOverall} ≔ \mathrm{DiffEquation}\left(\left[L\cdot i\'\left(t\right)\+R\cdot i\left(t\right) \= V\left(t\right)-K\cdot \mathrm{\θ}\'\left(t\right)\, J\cdot \mathrm{\θ}\'\'\left(t\right)\+b\cdot \mathrm{\θ}\'\left(t\right) \= K\cdot i\left(t\right)\right]\, \left[V\left(t\right)\right]\, \left[\mathrm{\θ}\left(t\right)\right]\right)\: \\ \mathrm{PrintSystem}\left(\mathrm{sysOverall}\right) \end{gathered}$$

To obtain the transfer function model of our system, we convert the overall differential equation model ($\mathrm{sysOverall}\)$ into a transfer function model using the $\mathrm{TransferFunction}$ command.  $\mathrm{TFopen}$ is the open loop transfer function for the angular position $\mathrm{\Θ}\left(t\right)$ from the input $V\left(t\right)$ and is defined as:

$$\begin{gathered} \mathrm{TFopen} ≔\mathrm{TransferFunction}\left(\mathrm{sysOverall}\right)\: \\ \mathrm{PrintSystem}\left(\mathrm{TFopen}\right) \end{gathered}$$

The state space model of a system is defined as follows:

For the case of the DC motor, the states $\mathrm{x1}\(t\)$, $\mathrm{x2}\(t\)$ and $\mathrm{x3}\(t\)$ represent current, $i\left(t\right)$, angular position, $\mathrm{\Θ}\left(t\right)$, and angular velocity, $\mathrm{\Θ}\'\left(t\right)$, respectively.

To obtain the state space model of our system, we convert the overall differential equation model ($\mathrm{sysOverall}\)$ into a state space model using the $\mathrm{StateSpace}$ command.  $$$$

$$\begin{gathered} \mathrm{SSopen} ≔ \mathrm{StateSpace}\left(\mathrm{sysOverall}\right)\: \\ \mathrm{PrintSystem}\left(\mathrm{SSopen}\right) \end{gathered}$$

The $\mathrm{ResponsePlot}$command can be used to plot the angular position $\mathrm{\Theta }\left(t\right)$with parameter set, $\mathrm{params}$, initial conditions, $\mathrm{ic}\,$ and input voltage signal, $\mathrm{Vin}\,$defined below:

$\mathrm{params}$

$\mathrm{ic}≔ \left[\mathrm{\θ}\left(0\right) \= 0\, \left(\mathrm{D}\left(\mathrm{\θ}\right)\right)\left(0\right) \= 0\, i\left(0\right) \= 0\right]$

$\mathrm{ResponsePlot}\left(\mathrm{TFopen}\,\mathrm{Heaviside}\left(t\right)\, \mathrm{output}\=\mathrm{\θ}\left(t\right)\, \mathrm{initialconditions}\=\mathrm{ic}\, \mathrm{parameters}\=\mathrm{params}\, \mathrm{duration}\=3\, \mathrm{thickness}\=2\,\mathrm{gridlines}\=\mathrm{true}\,\mathrm{title}\=''Change\; of\; shaft\; position:\; theta(t)''\, \mathrm{labels}\=\left[''time[sec]''\, ''theta(t)[rad]''\right]\, \mathrm{labeldirections}\=\left[\mathrm{horizontal}\,\mathrm{vertical}\right]\, \mathrm{font}\=\left[\mathrm{HELVETICA}\,10\right]\, \mathrm{labelfont}\=\left[\mathrm{HELVETICA}\,10\right]\right)$

The system can also be analyzed in the frequency domain. Substituting the parameter values into the open loop transfer function in (2.2.1) results in:$$

$\mathrm{eval}\left(\mathrm{TFopen}:-\mathrm{tf}\, \mathrm{params}\right)$

The above transfer function can be plotted as a magnitude plot and a phase plot.

$\mathrm{BodePlot}\left(\mathrm{TFopen}\, \mathrm{parameters}\=\mathrm{params}\, \mathrm{hertz}\=\mathrm{true}\, \mathrm{thickness}\=2\,\mathrm{labelfont}\=\left[\mathrm{HELVETICA}\,10\right]\, \mathrm{font}\=\left[\mathrm{HELVETICA}\,10\right]\right)$

Our objective is to be able to rotate the shaft position, $\mathrm{\Θ}\left(t\right)$, of the motor to specified desired positions.  In control terms, this translates to minimizing the motor shaft position error. As such, we put higher weighting on the position state, $\mathrm{x2}\(t\)$, in the Q matrix in our LQR controller design.

$$\begin{gathered} A≔\mathrm{eval}\left(\mathrm{SSopen}:-a\,\mathrm{params}\right)\; \\ B≔\mathrm{eval}\left(\mathrm{SSopen}:-b\,\mathrm{params}\right)\; \end{gathered}$$

$\mathrm{Qlqr}$

The controller gain for the full state feedback LQR controller is then computed using the custom procedure LQR( ). With the Q matrix defined in (3.1.2) and $\mathrm{Rlqr}\=0.1$the result is the gain vector $\mathrm{Kss}$ given by:

$\mathrm{Kss} ≔ \mathrm{LQR}\left(A\, B\, \mathrm{Qlqr}\, \mathrm{Rlqr}\right)$

The numeric solution can now be obtained for the step response of the open loop (uncontrolled) and closed loop (controlled) cases.

Open loop (uncontrolled) response Closed loop (controlled) response Input functions used examine the open loop and closed loop response are defined as: $\mathrm{InputSig}\;$ $\mathrm{PrintSystem}\left(\mathrm{DEopen}\right)$ $\left[\begin{array}{l}\mathbf{Diff.\; Equation}\\ \mathrm{continuous}\\ \mathrm{1\; output(s);\; 1\; input(s)}\\ \mathrm{inputvariable}\=\left[V\left(t\right)\right]\\ \mathrm{outputvariable}\=\left[\mathrm{\theta }\left(t\right)\right]\\ \mathrm{de}\=\{\begin{array}{l}\[\frac{\ⅆ}{\ⅆt}\mathrm{x1}\left(t\right)\=-\frac{R\mathrm{x1}\left(t\right)}{L}-\frac{K\mathrm{x3}\left(t\right)}{L}+\frac{V\left(t\right)}{L},\\ \frac{\ⅆ}{\ⅆt}\mathrm{x2}\left(t\right)\=\mathrm{x3}\left(t\right),\\ \frac{\ⅆ}{\ⅆt}\mathrm{x3}\left(t\right)\=\frac{K\mathrm{x1}\left(t\right)}{J}-\frac{b\mathrm{x3}\left(t\right)}{J},\\ \mathrm{\theta }\left(t\right)\=\mathrm{x2}\left(t\right)\]\end{array}\end{array}\right.$ (3.1.4) $$ The open loop response is defined as:   $\mathrm{PrintSystem}\left(\mathrm{sysOpenLoop}\right)$ $\left[\begin{array}{l}\mathbf{Diff.\; Equation}\\ \mathrm{continuous}\\ \mathrm{1\; output(s);\; 1\; input(s)}\\ \mathrm{inputvariable}\=\left[V\left(t\right)\right]\\ \mathrm{outputvariable}\=\left[\mathrm{\theta }\left(t\right)\right]\\ \mathrm{de}\=\{\begin{array}{l}\[\frac{\ⅆ}{\ⅆt}\mathrm{x1}\left(t\right)\=-\frac{R\mathrm{x1}\left(t\right)}{L}-\frac{K\mathrm{x3}\left(t\right)}{L}+\frac{V\left(t\right)}{L},\\ \frac{\ⅆ}{\ⅆt}\mathrm{x2}\left(t\right)\=\mathrm{x3}\left(t\right),\\ \frac{\ⅆ}{\ⅆt}\mathrm{x3}\left(t\right)\=\frac{K\mathrm{x1}\left(t\right)}{J}-\frac{b\mathrm{x3}\left(t\right)}{J},\\ \mathrm{\theta }\left(t\right)\=\mathrm{x2}\left(t\right)\]\end{array}\end{array}\right.$ (3.1.5) $$ The closed loop response is defined as: $\mathrm{PrintSystem}\left(\mathrm{sysClosedLoop}\right)$ $\left[\begin{array}{l}\mathbf{Diff.\; Equation}\\ \mathrm{continuous}\\ \mathrm{4\; output(s);\; 1\; input(s)}\\ \mathrm{inputvariable}\=\left[\mathrm{trim}\left(t\right)\right]\\ \mathrm{outputvariable}\=\left[\mathrm{x1}\left(t\right)\,\mathrm{x2}\left(t\right)\,\mathrm{x3}\left(t\right)\,V\left(t\right)\right]\\ \mathrm{de}\=\{\begin{array}{l}\[\frac{\ⅆ}{\ⅆt}\mathrm{x1}\left(t\right)\=-\frac{R\mathrm{x1}\left(t\right)}{L}-\frac{K\mathrm{x3}\left(t\right)}{L}+\frac{V\left(t\right)}{L},\\ \frac{\ⅆ}{\ⅆt}\mathrm{x2}\left(t\right)\=\mathrm{x3}\left(t\right),\\ \frac{\ⅆ}{\ⅆt}\mathrm{x3}\left(t\right)\=\frac{K\mathrm{x1}\left(t\right)}{J}-\frac{b\mathrm{x3}\left(t\right)}{J},\\ \mathrm{\theta }\left(t\right)\=\mathrm{x2}\left(t\right),\\ V\left(t\right)\=-0.523890948743337415\mathrm{x1}\left(t\right)-3.16227766016837997\mathrm{x2}\left(t\right)+3.16227766016837997\mathrm{trim}\left(t\right)-0.322243623661871481\mathrm{x3}\left(t\right)\]\end{array}\end{array}\right.$ (3.1.6) $$ The plot of the open loop state response is: $\mathrm{StatePlot1}$ $$ The plot of the closed loop state response is: $\mathrm{StatePlot2}$ $$ The plot of the open loop state response is: $\mathrm{OutputPlot1}$ The plot of the closed loop angular position $\mathrm{\Θ}\left(t\right)$ step response is: $\mathrm{OutputPlot2}$

In the controlled case, the shaft position of the DC motor follows the step change in the desired position (defined by the step input) and settles into the desired steady state value as expected.

We want to vary the weighting factor R for the LQR controller to see how it affects the system output response.

The R factor is the weighting coefficient for the input signal. By increasing R, we are putting more penalty on the control signal.  

As such, the resulting controller will try to minimize the control efforts used to achieve the desired output.

Consider the results obtained earlier for R = 0.1.

Next, we generate the controller and the corresponding step time response for each value of R from 0.1 to 2. To compare the performance of each of the controllers, we define the control effort of the controller as the square of the input signal, summed over time. This associated control effort (or energy-like) value for each controller is the total cost in operating the individual controllers.

The following plot illustrates the total cost (integral of the energy like control function) for different values of R.

$\mathrm{TotalCostPlot}$

To see the effect of varying R, you can drag the slider to select different values of R and observe the resulting step response and the control effort.