The rotary inverted pendulum control example is a case study to demonstrate modelbased design of cyberphysical systems. Controlling an inverted pendulum is a classic example to show how feedback control loop can be used to stabilize an inherently unstable system, and how parameterization of physical and cyber attributes can change the behavior of the system.
In this article we use the inverted pendulum example to go through the typical CPS design process step by step. These steps include mathematical modeling, control algorithm design, simulation software development and design space exploration using a simulationbased approach.
A specific research problem we can address with this example is managing power and energy in cyber physical systems. In this problem, the main question is how to select cyber system parameters to get the optimal power consumption. Similarly, we can study the effect of physical system parameters as well as cyber parameter to obtain a holistic model of the CPS containing notion of both actuation and processing power. In detail we discuss the following items:
 Inverted Pendulum Problem Definition
 Implementation
 Modeling
 Control Algorithm Design
 Simulation
 Design Space Exploration
 Related Publications
Inverted Pendulum Problem Definition
The primary goal of the inverted pendulum system is to stabilize the pendulum in an upright position. Therefore, a servo arm has to be moved appropriately to keep the pendulum angle close to zero with respect to vertical direction. Beside the stability we also consider two nonfunctional aspects in the problem definition: the control quality and the power consumption requirements.
We implemented the pendulum system in our lab and modeled the system with Simulink and Modelica.
Implementation
We implemented the inverted pendulum in our lab using the Batan S1213 R/C Servo and EKCLM3S6965 Texas Instruments ARM CortexM3 evaluation board. Pendulum angle and angular speed is measured by S1 Optical Shaft Encoder.
The source codes of the ARM microcontroller is available at "usebatanservo" branch of following Github repository:
hmirzai/invpend
Modeling
In order to design the control algorithm and also simulate the system, we need to model it first. Rotary inverted pendulum is simple enough to make us able to use direct equation derivation method. In this approach, the mathematical equations are derived using laws of Physics.
Rotary inverted pendulum consists of mechanical and electromechanical parts. Therefore, using laws of Dynamics and Electromagnetic we can obtain the required equations. First, the equations of subsystems, i.e. pendulum dynamics and servo motor are generated and then they will be put together to build the overall physical system model.
Dynamic Equations
We need to derive the differential equations which explains the relation between the variables and and also the effect of motor torque (the actuation variable) to the servo arm angle .
Using NewtonEuler method we can derive these equations as following *:
* To see the symbol definitions used in the above equations and rest of the page click on Symbol Definitions.


These equations are highly nonlinear and therefore hard to deal with using linear controllers. The good news is that we can work with the linearized version for the purpose of control algorithm design because when the system is stabilized, the variations of all of the variables in the equations are small and the nonlinear dynamic behavior is not unfolded.
Although the linearized dynamic equations should be used for the control algorithm design, the original nonlinear version can be used perfectly for the simulation purposes. To do this, we can solve the first two equations for variables and to get the statespace equations in the following form:
These equations can be implemented easily in simulation tools like Matlab Simulink or Modelica.
To see the details of dynamic equation derivation process click on Derivations.
Assuming the mass of pendulum rod is negligible, we can write the NewtonEuler equations for the rigid body consisting of pendulum weight and rod in the coordinate frame as following:
\begin{align}
\sum \mathbf{F} &= m \mathbf{a}_G \label{eq:nef} \\
\sum \mathbf{M}_A &= \mathbf{r}_{G/A} \times m \mathbf{a}_G + \frac{d}{dt} \mathbf{H}_A, \label{eq:nem}
\end{align}
Using the definition
\begin{align}
\mathbf{s} = \mathbf{r}_{G/A} \times m \mathbf{a}_A ,
\end{align}
\eqref{eq:nef} and \eqref{eq:nem} can be expanded for the three coordinate axes as following:
\begin{align}
\sum F_x &= m a_x, \\
\sum F_y &= m a_y, \\
\sum F_z &= m a_z, \\
\sum {M_A}_x &= I_{xx} \alpha_x  (I_{yy}  I_{zz}) \omega_y \omega_z + s_x , \label{eq:nemx}\\
\sum {M_A}_y &= I_{yy} \alpha_y  (I_{zz}  I_{xx}) \omega_z \omega_x + s_y ,\\
\sum {M_A}_z &= I_{zz} \alpha_z  (I_{xx}  I_{yy}) \omega_x \omega_y + s_z ,
\end{align}
in which are the components of inertia tensor matrix. The offdiagonal component of that matrix are all zero and the diagonal components are:
\begin{align}
I_{xx} = mL^2 , \quad I_{yy} = mL^2 , \quad I_{zz} = 0. \label{eq:inertiatensor}
\end{align}
Next, we will calculate the terms appeared in above equations in terms of angles and and their derivatives using kinematic equations. The auxiliary coordinate frames and are used for this purpose (see the Figure).
\begin{align}
\boldsymbol\omega &= \dot{\alpha} \mathbf{K} + \dot{\theta} \mathbf{i} \implies \boldsymbol\omega = \dot{\theta} \mathbf{i} + \dot{\alpha}\sin \theta \mathbf{j} + \dot{\alpha}\cos \theta \mathbf{k} \\
\omega_x &= \dot{\theta},\quad \omega_y = \dot{\alpha}\sin \theta,\quad \omega_z = \dot{\alpha}\cos \theta \\
\alpha_x &= \ddot{\theta},\quad \alpha_y = \dot{\theta}\dot{\alpha}\cos \theta + \ddot{\alpha} \sin \theta, \alpha_z = \dot{\theta}\dot{\alpha}\sin \theta + \ddot{\alpha} \cos \theta \\
\mathbf{s} &= L\mathbf{k} \times m(l\dot{\alpha}^2\mathbf{i}+l\ddot{\alpha}\mathbf{j'}) = mLl\dot{\alpha}^2 \mathbf{j'}  mLl\ddot{\alpha}\cos \theta \mathbf {i} \label{eq:vec_s}
\end{align}
Substituting the terms in \eqref{eq:nemx} by what we have in \eqref{eq:inertiatensor} to \eqref{eq:vec_s} and adding up all the momentum applied by different forces on the pendulum rigid body system we will get
\begin{align}
\mu_p \dot{\theta} + mgL\sin \theta = mL^2 \ddot{\theta}  mL^2 \dot{\alpha}^2 \sin \theta \cos \theta  mLl\ddot{\alpha}\cos \theta \label{eq:thdd}
\end{align}
Since we have two variables, we need another equation in addition to \eqref{eq:thdd} to build the statespace (SS) model. Also, we need to involve the torque applied by the actuator in the new equation. For this purpose, torque equation around the axis of the can be written as following
\begin{align}
J \ddot{\alpha} = \tau  \mathbf{K}.(l\mathbf{i}\times\mathbf{F})  \mu_m \dot{\alpha}. \label{eq:aldd}
\end{align}
Then, we need to calculate using \eqref{eq:nef} and substitute the result in \eqref{eq:aldd}. To do this, first we need to obtain using following equation:
\begin{align}
\mathbf{a}_G =& \mathbf{a}_A + \alpha \times \mathbf{r}_{G/A} + \boldsymbol\omega \times (\boldsymbol\omega \times \mathbf{r}_{G/A}) \nonumber \\
=& (l \dot{\alpha}^2 + L \sin \theta \ddot{\alpha} ) \mathbf {i}
+ (l\cos\theta \ddot{\alpha}  L \ddot{\theta} + L \sin \theta \cos \theta \dot{\alpha}^2) \mathbf{j} \nonumber \\
&+ (l\sin \theta \ddot{\alpha}  L \sin^2 \theta \dot{\alpha}^2  L \dot{\theta}^2) \mathbf{k} \nonumber \\
=& F_x \mathbf{i} +F_y \mathbf{j} + F_z \mathbf{k} \label{eq:F}
\end{align}
Now we can use \eqref{eq:F} to rewrite \eqref{eq:aldd} as following
\begin{align}
J \ddot{\alpha} = \tau  lF_y \cos \theta + lF_x \sin \theta  \mu_m \dot{\alpha}. \label{eq:aldd2}
\end{align}
Actuator Modeling
In our example, the Actuator is a geared DC motor. Following diagram shows the widely used for DC motor which is used to derive the actuator equations
The equations that describes the actuator dynamics are
\begin{align*}
L_a \dot{i}_a + R_a i_a &= v_m  K_b \dot{\alpha}_m \\
\tau_m &= K_t i_a  B_m \dot{\alpha}_m
\end{align*}
Also, assuming the gearbox attached to the DC motor is ideal we will have
\begin{align*}
\tau &= \tau_m n \\
\alpha &= \alpha_m / n
\end{align*}
where is the gear ration.
These above equations are linear and we neglect the nonlinearities of DC motor dynamics.
Control Algorithm Design
With the models of physical plant and actuator we can design the controller. The overall system model including the controller is like the block diagram shown here. All the physical components including the actuator (dc motor and gearbox) and pendulum dynamics are shown inside the dashed box. The controller sees the whole thing as single component.
LQR controller has been used for stabilizing the inverted pendulum. Actually, it is also used in many similar works. The main objective of LQR design is to minimize a linear combination of control performance and an energy consumption metrics. The coefficients can be tuned to get the required tradeoff between those metrics. One example of a const function is
\begin{align*}
J = \int_0^\infty \left( q\alpha\alpha_r^2 +v_m^2 \right) dt \label{eq:cost2}
\end{align*}
represents the tracking error , that is the difference between the actual servo arm angle and the commanded angle. Therefore, it is a metric of control performance. is the squared value of the motor terminal voltage which is roughly proportional to the motor power consumption. The constant q can be used to adjust the tradeoff between the control performance and power consumption.
From LQR theory we know that the control law will be the simple proportional state feedback. In our example it means that we should apply the motor voltage using the following equation:
\begin{align*}
v_m =  K \mathbf{x},
\end{align*}
is the vector containing all the system states and is a matrix gain. can be calculated using offtheshelf packages like lqr function of Matlab Control toolbox. Before doing that we should linearize the system and express in the standard linear statespace form. We also need to combined the statspace models of components inside the dashed block (actuator and pendulum dynamics into a single statespace model.
To see the linear statespace equation of our example click on linear statespace model.
Linearized StateSpace model of pendulum dynamics:
\begin{align}
\dot{\mathbf{x}}_p = A_p \mathbf {x}_p + B_p u_p, \label{eq:ss}
\end{align}
where
\begin{align}
\mathbf{x}_p &=(\begin{array}{cccc}\alpha & \dot{\alpha} & \theta & \dot{\theta}\end{array})^T, \quad u_p = \tau \\
A_p &= \left. \frac{\partial \mathbf{f}}{\partial \mathbf{x}_p} \right_{\mathbf{x}_p=\mathbf{x_{p0}}=\mathbf{0}} = \left( \begin{array}{cccc}
0 & 1& 0& 0 \\
0 & \frac{\mu_m}{J}& \frac{glm}{J}& \frac{l\mu_p}{JL} \\
0 & 0& 0& 1 \\
0 & \frac{l\mu_m}{JL}& \frac{g(ml^2 + J)}{JL}& \frac{\mu_p(ml^2 + J)}{JL^2m}
\end{array} \right) \\
B_p &= \left. \frac{\partial \mathbf{f}}{\partial u_p} \right_{u_p=u_{p0}=0} =
\left( \begin{array}{cccc} 0 & \frac{1}{J} & 0 & \frac{l}{JL} \end{array} \right)^T
\end{align}
Linearized StateSpace model of actuator (DC motor):
\begin{align}
\dot{x}_a &= A_a x_a + B_a \mathbf{u}_a = \frac{R_a}{L_a} x_a + (\begin{array}{cc} \frac{1}{L} & \frac{K_a}{L_a} \end{array}) \mathbf{u}_a \label{eq:motss1} \\
y &= C_a x_a + D_a \mathbf{u}_a = K_t x_a + (\begin{array}{cc} 0 & B_m \end{array}) \mathbf{u}_a \label{eq:motss2}
\end{align}
where
\begin{align}
x_a &= i_a,\quad\mathbf{u}_a = (\begin{array}{cc}v_m & \dot{\alpha}_m\end{array})^T,\quad y = \tau_m \\
\end{align}
Simulation
Pendulum Angle(deg) vs Time(s)
Servo Angle(deg) vs Time(s)
3D animation shown above demonstrates the results of simulation. The simulation model is built using both Simulink and OpenModelica. Click on "Free Fall" and "LQR Control" button to see the two different scenarios:
 Free fall: in this scenario, no torque is applied by the actuator to the pendulum. Therefore, just a small deviation from vertical angle causes the pendulum to fall down. Thanks to nonlinear equations we can observe the large signal behaviour of the physical system.
 LQR Control: in this scenario, the LQR controller is in place and keeps the pendulum in upright position while moves it from one place to another back and forth. The square wave signal in the servo angle plots shows the commanded signal and the yellow curve shows the actual value which is the result of closedloop control. Although the angular deviations are relatively small, the nonlinear simulation is still useful because it can be used to verify the goodness of approximation caused by linearizing the original nonlinear equations.
We can use the simulation not only to verify the design process, but also for design space exploration and fine tuning the system design done so far or to explore the dimensions which are not included in the initial system design. For example, all the design process at previous stages include only continuoustime modeling of the whole system (physical system + controller). But, in CPS we are implementing the continuoustime controller using a discretetime cyber system. Simulation model gives us the opportunity to study all the imperfections caused by discretization process.
Design Space Exploration
The simulation software developed in previous stages can be used as a tool to explore the design space. For this reason, we can simply run the simulation software for every setting point in the design space. However, since the design space grows exponentially as we add more dimensions to the problem, running the simulation model is costly. To reduce the computation time, we can consider two solutions to mitigate this issue. First, we should try to produce the fastest possible version of the simulation software. Second, we should leverage parallel computing tools because different setting points in the design space are independent of each other.
We have used OpenModelica to develop simulation software and the increase the speed of each simulation run. OpenModelica translate the simulation model to native C++ software which causes considerable speed up as opposed to interpreterbased solutions like Matlab Simulink. Although Simulink also comes with a compiler that can translate the simulation model, our tests showed that Modelica provides better performance and flexibility especially when we change the sampling time parameter of the controller subsystem.
To realize the second solution, i. e. to get effective use of parallel computing platforms, we have used multiprocessing package in the Python script which is developed to run the simulation software in nested iterations of each design dimension. By doing this, we can make use of all the CPU cores available on the parallel processing platform.
The explained methodology can be used to answer many CPS design question. Two examples of these questions are explained as follows:
What is the relation of the pendulum length to the controller clock speed?
From experience with pendulums with different length we know intuitively that longer pendulum length results in slower physical dynamics. Therefore we expect that we can increase control task period which leads to lower CPU frequencies. But, what is the exact relation.
To get the answer to this question we can run the experiment for different values of length and control task period and see how the system response changes. The figure shown above shows the result of such analysis. For a fixed value of pendulum length, the system goes toward instability as we increase the control task period. Same thing happens for a constant control task period as we decrease the pendulum length. The region of design space in which we have relative stability is shown roughly in the figure.
What is the optimal setting for the control task period to get the minimum energy consumption?
The overall energy consumption is the sum of actuation energy and processing energy. We saw that the control design cost function somehow considers the actuation energy but for two main reasons we need to extend our model to support following properties of the actual system: 1 The cost function represents the continuoustime actuation energy usage, but after discretizing the system (running the controller as a embedded system task with some period) it is not a valid energy expression. 2 The dcmotor power is not exactly proportional to the squared value of terminal voltage because it is not a pure resistive load.
We can use the simulation to calculate the energy using more accurate equation comparing to what is considered the LQR cost function. Also, we can add the processing energy to the total energy consumption in the simulation model. For shorter control task periods we expect better control performance and less actuation energy consumption. In the other hand processing power is higher for shorter task periods because we have to invoke the control procedure more frequently in unit of time in that case. This intuition suggests that there should be an optimal setting for control task period that minimizes the total energy consumption. Following figure made by running the simulation software for different values of the control period (sampling time) shows that this intuition is true in our example. In this figure, is the weight parameter in LQR controller cost function. (see Control Algorithm Design in this page)