### 2. The Seesaw/Pendulum Process

In this section, we will illustrate how the symbolic capabilities of Mathematica can be used to derive a model of a nonlinear system. Using MathCode C++, the large symbolic expressions in the nonlinear model are then converted to C++ code and used to simulate the system very efficiently. The system is a laboratory process frequently used in control education known as the seesaw/pendulum process. One implementation of the seesaw/pendulum process has been developed by Quanser Consulting (www.quanser.com).

#### The System

The process consists of a seesaw, two carts called and , two parallel tracks, an inverted pendulum, and a weight. Each cart can be driven by a DC motor controlled by an input voltage. Cart carries the weight and cart carries the inverted pendulum attached by a friction-free joint. The carts can be moved along the tracks by controlling the input voltages of the DC motors.

We start by modeling the open loop system, that is, without any feedback from measured signals. The forces and acting on each cart are chosen as inputs. We will use the Lagrangian methodology to obtain a nonlinear model of the system and then linearize it. The linearized model is used for control system design where the Control System Professional application package is used. The closed-loop system is then simulated both within Mathematica and by external code generated using MathCode C++.

Figure 1. The pendulum/seesaw process.

In Figure 1 we introduce names of variables and constants describing the pendulum/seesaw process.

Variables:
translation of cart 1 from center of track
translation of cart 2 from center of track
angle of seesaw with vertical
angle of inverted pendulum with normal to track
force applied to cart 1
force applied to cart 2

Constants:
J inertia of seesaw with track at height h
mass of seesaw with track
c height of center of gravity of the seesaw from pivot point
h height of track from pivot point
mass of cart 1 (weight cart, on the back track)
mass of cart 2 (pendulum cart, on the front track)
mass of pendulum
center of mass of pendulum rod (half of full length)
g gravitational acceleration

The physical values of the above constants are stored as rules in Mathematica to be used in later simulations.

#### Modeling

Following the Lagrangian methodology [5], the system is divided into a number of subsystems whose potential energy and kinetic energy are computed in terms of the generalized coordinates introduced in Figure 1. The Lagrangian, which is the difference between the total kinetic and potential energy, can then be used to derive the equations of motion for the system.

Computation of the Lagrangian

Coordinates of the center of mass of the seesaw:

The potential and kinetic energies of the seesaw:

Coordinates of the center of track:

Coordinates of cart 1:

The potential and kinetic energies of cart 1:

Coordinates of cart 2:

The potential and kinetic energies of cart 2:

Coordinates of the center of mass of the pendulum:

The potential and kinetic energies of the pendulum:

The total potential energy:

The total kinetic energy:

The Lagrangian of the system:

becomes

Computation of the Equations of Motion

The equations of motion are given by the partial differential equations that must satisfy

where is the Lagrangian, are the generalized coordinates, and are the generalized forces associated with each coordinate. In this case, the quadruple corresponds to and , , . We derive each equation and simplify it.

These equations are coupled ordinary differential equations (ODEs) of second-order in the generalized coordinates . To rewrite these in standard state-space form (a system of first-order ODEs) we introduce the following states:

The inputs to the system are

We observe that the second-order time derivatives of the positions and angles appear linearly, which makes it easy to solve for these in terms of the states. Second-order time derivatives will always appear linearly in the equations derived from the partial differential equation (6) that the Lagrangian of the system has to satisfy. This is due to the fact that the Lagrangian only consists of at most first-order time derivatives and the second-order derivatives appear according to the chain rule when differentiating the term in equation (6) with respect to time.

Solve for second-order derivatives.

These solutions will become a part of the right-hand sides of the differential equations used for simulating the system.

We introduce conversion rules to get rid of explicit time dependence. This makes some expressions look simpler and takes less space.

We compute the right-hand side of the nonlinear state-space form (equation (1)) of the equation

which becomes very large! Observe that the first four entries in f correspond to pure integrations. This allow us to write the model as a system of first-order differential equations as desired.

The Linearized Model

Since we will design a controller using methods for linear systems, we need to linearize the nonlinear state space model of the system around the origin.

This can easily be obtained using the `Linearize` function in Control System Professional.

We extract the matrices from the linearization and substitute with the physical values of the parameters to get numerical entities.

The linear model of the system allows us to use standard methods for controller design.

#### Controller Design

The method for computing the feedback law requires the linear model to be controllable. This implies that the system can be stabilized and that the method will find a numerical instantiation of such that this is obtained. To check controllability we use the `Controllable` command in Control System Professional.

Compute a full state feedback controller using the LQG-method (`LQRegulatorGains` is a Control System Professional command).

The control law to be used is , where are measurements of the states. This gives the following closed loop system to simulate .

#### Simulation and Code Generation

Preliminaries

We store the states as a vector.

We close the feedback loop, which gives the following right-hand side of the state-space equations.

Change the names of variables to only ASCII characters (needed for code generation).

The right-hand side of the closed-loop state-space equations (intended for code generation).

We define the function `ClosedLoopPendulumEquationsRHS`, which has the huge right-hand side expression from above as its body. The type declarations `Real` and `-> Real[8]` below are the syntax in MathCode C++ for providing type information. Explicit type declarations are needed to generate efficient C++ code.

Simulation of the Nonlinear Model within Mathematica

We define to be a short notation for the `ClosedLoopPendulumEquationsRHS`.

The state vector is defined.

The differential equations for the closed-loop system is given.

We need to specify some initial conditions for the simulation. Assume that cart 1 is positioned at -1, the seesaw is horizontal, cart 2 is positioned at +1, and the pendulum has a small deviation of 5.7° from the vertical axis. Furthermore, assume that the system is at rest at these positions (all time derivatives are zero).

Nonlinear System Simulation

We simulate the system using `NDSolve`.

Linear System Simulation

We simulate the linear system using `NDSolve`.

Plot the result for the first four state variables from the nonlinear simulation together with the result for the linear system.

Figure 2. Initial value response for the linear (red) and nonlinear (blue) systems.

We notice that the linear response differs significantly from the nonlinear one. The difference between the linear and nonlinear response will be much smaller if the initial values are chosen to be smaller since the linear model only is a good approximation for small values of the states and inputs.

External Simulation of the Nonlinear Model

To generate simulation code to be run outside Mathematica, we have to provide a differential equation solver that can be compiled or linked when the executable program is generated. In this case, we have chosen to implement a very simple Runge-Kutta solver in Mathematica, which will be included in the generated code. Another solution could be to use a solver from a publicly available software library like the solvers at www.netlib.org.

Code for External Simulation of the Nonlinear Model

Here follows an implementation of a state equation solver using the Runge-Kutta method.

We use `RK` in the following ODE solver.

To get a correct estimate of the time spent by the external code for simulating the system, we define the following functions that repeat the external simulation n times.

Compilation

We are now ready to compile the package. The right-hand side of the state equations are optimized using common subexpression elimination.

The corresponding code can be inspected in a text editor.

All necessary files for building binaries are now generated. The `MakeBinary` command can build either a standalone application using the code specified by an option `GenerateMainFileAndFunction` or a MathLink version, which can be easily installed into Mathematica. We build the MathLink version.

To use the generated code for the simulation, we only have to install the code using the `InstallCode` command from MathCode C++.

External Simulation: Uncompiled Code (Computations within Mathematica)

We plot the result to verify the Runge-Kutta solver.

Figure 3. Initial value response for the nonlinear system computed internally by Mathematica.

External Simulation: Compiled Code (External Computations)

We install the executable code.

The external simulation is performed 200 times to get an accurate measure of its timing.

We plot the result to verify the compiled Runge-Kutta solver.

Figure 4. Initial value response for the nonlinear system computed by the external code.

The responses in Figure 3 and Figure 4 are identical.

Performance Comparisons

The difference in performance between the uncompiled and compiled version of the ODE solver is 1419.03.

A somewhat unfair comparison between the timings of the built-in ODE solver and the compiled solver gives 15.8667.

#### Clean-up

Uninstall the code and delete the temporary files.

Remove the functions for which we have generated code.

Unset the symbol f storing the right-hand side of the state-space equations.