Jtam-A4.dvi


JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

51, 4, pp. 1027-1038, Warsaw 2013

CONTROL SYSTEM DESIGN FOR MULTI-ROTOR MAV

Wei Wang, Hao Ma

College of Information and Control, Nanjing University of Information Science and Technology, Nanjing, China

Chang-Yin Sun

School of Automation, Southeast University, Nanjing, China

w-mail: wwcb@nuist.edu.cn

Abstract MAVs (Micro Air Vehicles) have attracted more and more attention in recent
years because of their wide applications, such as investigations of inhospitable environment,
nuclear or hazardous material polluted regions and other inaccessible areas. In this paper,
both attitude and altitude controls of a Quad-Rotor typeMAVare discussed and analyzed.
To stabilize the attitude, one order memory based controller is selected, the performance
of which is independent of the detail model information, leading to enhanced robustness.
While for the altitude, aLinearQuadraticGaussian (LQG)controller is proposed.The expe-
rimental results of the real flight demonstrate the stability and robustness of the controllers.

Key words: MAV, attitude control, altitude control, LQG, memory based control

1. Introduction

Recently, a new class of flight vehicles, called Micro Air Vehicles (MAVs) which have great
potential for military, industrial as well as civilian applications, has gained tremendous atten-
tion. Since the 1990s, with the development of the Micro Electro Mechanical System (MEMS)
technology, a single rotor helicopter has been the favorable platform for MAV development.
However, the disadvantages such as complexity, instability and limited payload slowed down
its development. While the multi-rotor vehicle, especially the Quad-Rotor vehicle, is studied
in more research experiments. The Quad-Rotor MAV, compared with a single rotor helicopter
and dual-rotor vehicle, has excellent capabilities of small size, superior payload, low cost and
simple manufacture. As a result, it can be widely used in military reconnaissance, meteorolo-
gical survey, environmental monitoring and air transportation in emergency and other aspects
(Martinez, 2007).
Normally, MAVs need to fly autonomously to accomplish those difficult missions. However,

many technical difficulties are induced by their small sizes and structures, especially the in-
stallation of communication apparatus and sensors. Until now, tremendous achievements have
beenmade in solving these difficulties bymany groups and institutions such asMIT,Ascending
Technology, Microdrones, Draganfly, Heudiasyc University, Pennsylvania University and Chiba
University (Achtelik etal, 2009; Bachrach etal, 2009; [2], [5],Wangetal, 2008;Nonami etal, 2010).
To realize an autonomous flight, various types of control techniques have been proposed, inc-
luding back stepping (BS), nonlinear control, robust control and so on (Kendoul, 2012, Tayebi
and McGilcary, 2006; Xia etal, 2012). Recently, the discussion of flight control has been driven
to amature stage and the focus ofMAV research has been transferred to the vision based flight,
indoor obstacle avoidance (Hoffman etal, 2011), SLAM based control (Rudinetal, 2011) or the
development of new structure MAVs. Despite that, the performance of all flight controllers is
limited by the external disturbance and model errors. Therefore, only in certain circumstances
can theywell operate.Maybe it is not serious for research study, butwhen it is used for practical



1028 W.Wang et al.

application, the environmental restrictionwill be a remarkable imperfection. So it is still of great
significance to develop a really effective robust controller.
In this paper, themain considerations of systemdynamics of attitude andaltitudemotions as

well as the control algorithms are discussed.With the control input and output of the real plant
obtainedbyflight experiments, themathematicalmodel canbeobtainedby system identification
tools inMatlab.Toachieve theattitude control,weapplyamemorybasedcontrollerwhichmakes
use of certain historical information to build the control scheme. The greatest attractiveness of
this algorithm is that the control performance is not related to the detail of system dynamics
and outer disturbances, so that significant robustness is maintained. As to the altitude, we use
a barometer and accelerometer instead of the widely used GPS to improve the precision. Then
a model based optimal controller is applied. A Kalman filter is also applied to filter the noise
and estimate the unobserved state variables.

2. Controlled object

Our research platform is a 4-channel radio controlledQuad-RotorMAVdesigned byour research
group. The vehicle has a carbon fiber airframe, which provides great stiffness and rigidity while
keeps the weight down. The lift results from fourAKEbrushlessmotors. In order to drive these
brushless motors, Pentium-18A Ultra-PWM brushless driver which can support 50Hz-500Hz
PWM signal communication is selected. With the brushless driver, the maximum speed can
reach about 8000 rpm and payload of about 1kg can be acquired at the same time with four
propellers of APC1245. The flight time is about 10 minutes with a three-cell-2100mAh LiPo
battery without payload. Figure 1 shows an overview of our platform. And the specification is
shown in Table 1.

Fig. 1. Quad-Rotor typeMAV

Table 1. Specification of Quad-Rotor typeMAV

Parameter Value Unit

Maximum diameter 500 mm

Height 250 mm

Total mass 1.1 kg

Maximum lift 1.0 kg

The embedded control systemconsists of one IMUsensorXSENS, one triaxial accelerometer,
one barometer, one 12 bit ADC chipMCP3204 and one Arm7micro controller. The navigation



Control system design for multi-rotor MAV 1029

data for attitude control are read from theXSENS through anRS-232 serial communication. To
improve the accuracy of altitude information, a triaxial accelerometer ADXL335 from Analog
Devices Corporation together with an MPX5100Ap barometer is supplied to the Kalman filter
for estimation. The Arm7 receives control instruction by a 72MHz PPM receiver and sensor
data to run independent control loops for each axis. The entire processing is donewith a control
loop frequency of 400Hz. The configuration of the overall system can be seen in Fig. 2. In order
to obtain the flight data and prepare for the future research, communication between the vehicle
and the ground control station is accomplished at the frequency of 50Hz via a wireless modem
XBee-pro which can communicate up to 1.6kmwith its maximum power of 63mW.

Fig. 2. Configuration of Quad-Rotor system

3. Attitude modeling and controller design

This section mainly illustrates the attitude modeling and controller design. Generally, for the
X and Y directions, we divide the autonomous flight control loop into three loops: attitude,
velocity, and position. Therefore, the first step is to accomplish the attitude control, which
directly affects the flight performance.

In viewof the inevitable nonlinearity and the uncertain outer disturbances, it is impossible to
generate the accurate attitudemodel. Therefore, it is of great importance to design a controller
with strong robustness. In this section, a one order memory based controller is recommended.
The fundamental idea of the control algorithm is to utilize certain memorized information such
as the current system response, previous system response and past control experience to build
control action. The main advantage of this approach is that the controller performance has no
relation to the detail of system dynamics. In addition, the overall required memory space does
not grow with time and is much smaller thanmost existing learning control methods.

3.1. Attitude modeling

Tuning of rotational speed of every rotor leads to every degree of freedom of the Quad-
-Rotor, therefore, the coupled motion is very obvious. However, it is quite difficult to analyze
such a multiple-input andmultiple-output (MIMO) system. For simplification, the modeling of
the nonlinear attitude motion is carried out on the premise of neglecting the coupling in the
following discussion.



1030 W.Wang et al.

Plus, the attitude model of Roll is regarded as the same as that of Pitch because of the
symmetry of theQuad-Rotor structure. Thus, the attitudemodeling is conducted aiming at the
Roll angle only.
First, assume that the Quad-Rotor is a rigid body evolving in a three-dimensional space.

In the body coordinate B, the equations of motion for a vehicle of mass m subject to external
force F∈ R, torque τ ∈ R can be given as (Kendoul et al., 2007, 2008)

mV̇B +Ω×mVB =F JΩ̇+Ω×JΩ= τ (3.1)

where V = [u,v,w] and Ω = [p,q,r] represent the linear and angular velocities in the body-
fixed reference frame. The translational force F is a combination of gravity, main thrust and lift
forces from rotors.
Using Euler angles η = [φ,θ,ψ] and the convention ZY X, the rotation matrix from B to

the inertial reference frame I is expressed as follows

R=







cosθcosψ sinφsinθcosψ− cosφsinψ cosφsinθcosψ+sinφsinψ
cosθ sinψ sinφsinθsinψ+cosφcosψ cosφsinθsinψ− sinφcosψ
−sinθ sinφcosθ cosφcosθ






(3.2)

Note that the matrix R satisfies R−1 = RT and det(R) = 1. And considering the trans-
formation between the two related coordinates as seen in Fig. 3, it is possible to separate the
gravity from other forces and write the translational dynamics in I as follows

ξ̇=v mv̇=RFb −mgei3 (3.3)

with ξ = [x,y,z] and v = [ẋ, ẏ, ż) are the vehicle position and velocity in I. In the above
equation, g is the gravitational acceleration and Fb is the resulting force (without gravity)
acting on the airframe.

Fig. 3. Rigid body dynamics and the associated frames

Up to now, different techniques have been proposed to express the rotational dynamics
such as the quaternion, rotation matrix R or Euler angles. In the following, the Euler angles
parametrization is employed to express the rotational dynamics in an appropriate form for
control design.
To begin with, let us consider the kinematic relation between Ω and η̇

η̇=Φ(η)Ω (3.4)

with the Euler matrix

Φ(η)=







1 sinφtanθ cosφtanθ
0 cosφ −sinθ
0 sinφsecθ cosφsecθ






(3.5)



Control system design for multi-rotor MAV 1031

Note that the matrix Φ has singularities at θ =±π/2, and its inverse matrix Ψ(η)=Φ−1(η)
is given

Ψ(η)=







1 0 −sinθ
0 cosφ cosθsinφ
0 −sinφ cosθsecφ






(3.6)

BydifferentiatingEq. (3.4)with respect to time, the secondequationofEq. (3.1) canberewritten
as

η̈= Φ̇Ω+ΦΩ̇= Φ̇Ψη̇−ΦJ−1sk(Ω)JΩ+ΦJ−1τ (3.7)

The sk operation is defined here from R3 to R3×3 such that sk(x) is a skew symmetricmatrix
associated with the vector product sk(x)y := x×y.
As an assumption, the attitude angle is small enough to ignore the coupling, Φ and Ψ

become identity matrices. Therefore, the differential equation of Roll motion can be obtained

φ̈ = fφ(φ̇,φ)+J
−1
φ
τφ (3.8)

with Jφ and τφ are the rotational inertia and torque around x axis, respectively, and fφ(φ̇,φ)
represents a function about φ̇ and φ, which can be generated from Eq. (3.7).
Generally, the control input u to the actually obtained torque τφ can be thought as the first

order delay system, which can be thought as a nonlinearity or disturbance to be absorbed to
fφ(φ̇,φ). Then, we arrive at

φ̈ = f(φ̇,φ)+J−1
φ
u (3.9)

In thedifferential equationabove, theunknownparameter Jφ canbedeterminedbyexperiments,
while the detail information of f(φ̇,φ) still keeps unavailable. However, the control algorithm
to be discussed does not involve the specific expression of f, so it is out of necessity to conduct
a deep study about this function.

3.2. Memory based controller design

For attitude controller design, considering the nonlinearity and the strong coupling of real
Quad-Rotor motion, it is of great importance to design the controller with good robustness in
both stability and performance. Therefore, the memory based control technique is adopted to
enhance the robustness (Song, 2006).
First, rewrite the differential equation

φ̈ = f(φ̇,φ)+J−1
φ
u (3.10)

Defining the error e = φ−φ∗, the error system of the Roll model can be described as

ë = f(φ̇,φ)+J−1
φ
u− φ̈∗ (3.11)

where φ∗ is the reference input angle.
For the convenience of processing, a new variable s is introduced

s =(q +k0)e where k0 > 0 q =
d

dt
(3.12)

The advantage of this treatment is that the original 2-order tracking problem in respect of φ
is reduced to 1-order stabilization in s. Even if the system equation grows to n-order, this
processing approach can also be applied by just adjusting s to be s =(q+k0)

n−1e.



1032 W.Wang et al.

In order to substitute ë, we only need to differentiate s once for the input, leading to

ṡ =(q +k0)ė = ë+k0ė (3.13)

Replacing ë above by Eq. (3.11), we get

ṡ = f(φ̇,φ)− φ̈∗+k0ė+J
−1
φ
u (3.14)

Now, define the control input as

J−1
φ
u = v+η(·) (3.15)

where η(·) = φ̈∗ −k0ė, and v is a memory-based compensator accounting for system nonline-
arities and other disturbances. The closed-loop system in terms of s becomes

ṡ = f(φ̈,φ)+v (3.16)

By Euler differentiation, we get

sk+1 = sk +(fk +vk)T (3.17)

Here, T is the sampling period. Now performing backward time-shift in the above equation, we
have

sk = sk−1+(fk−1+vk−1)T (3.18)

Subtracting Eq. (3.18) from Eq. (3.17), gives

sk+1 =2sk −sk−1+(fk +vk)T − (fk−1+vk−1)T (3.19)

Now define vk as

vk = w0vk−1+w1sk +w2sk−1 (3.20)

with w0 =1, w1 =−2/T , w2 =1/T . Then, Eq. (3.19) can be rewritten as

sk+1 =(fk −fk−1)T = T∆fk (3.21)

For a fairly small T , ∆fk ≈ T(df/dt) |t = kT , so we have

|sk+1| ¬ T
2λ (3.22)

λ in the above equation is defined as

λ =max
(
∣

∣

∣

df

dt

∣

∣

∣ |t =0,T,2T,. . . ,nT
)

(3.23)

Considering the definition of s which can be regarded as a sliding hyperplanewith a chattering
range of T2λ, when s is convergent, the tracking error e and its derivative ė are also convergent.
The equations above indicate that if the variation of function f is not infinitely fast, the

tracking error and its derivative would be considerable small (since T is very small in general).
To summarize, for the above attitude model with a sampling period of T and reference

input φ∗, the memory based controller is designed as (by subtracting uk−1 = Jψ(vk−1+ηk−1)
from Eq. (3.15))

uk = uk−1−
2Jφ
T
(ėk +k0ek)+

Jφ
T
(ėk−1+k0ek−1)+Jφ(ηk −ηk−1) (3.24)

with Jφ being the rotational torque, η(·) = φ̈
∗−k0ė.



Control system design for multi-rotor MAV 1033

For simplification, η(·) can bemerged into f(·), so the control input expression is driven to

uk = uk−1−
2Jφ
T
(ėk +k0ek)+

Jφ
T
(ėk−1+k0ek−1) (3.25)

To verify the effectiveness of the controller, we carried out a series of simulations. For the
attitude system of Eq. (3.10), the detailed information of function f(·) keeps unknown, so it is
quite necessary to change f(·) in different forms

f1(φ̇,φ)= φ̇+φ+1 f2(φ̇,φ)= φ̇φ+ φ̇+φ+1 (3.26)

The simulation results of the controller are shown in Fig. 4. We can see that for different form
of function f(·) the memory based controller can keep excellent tracking performance, so that
the good robustness is realized.

Fig. 4. Simulation results of memory based controller (T =0.02s, K0 =1, J =0.173)

4. Altitude modeling and controller design

The altitude results from the combination force of four rotors, so the model can be easily
obtained. Then the LQG controller based on optimal controller is designed.

4.1. Altitude modeling

The altitude of the vehicle is controlled by the rotational speed of four rotors. According to
the previous work on the rotor performance test, the relationship between the control input u
and generated thrust F can be supposed as one order system

F(s)=
k

1+T1s
u(s) (4.1)

Here k and T1 are the constant parameters of the one order system. The relations between
thrust F and position can be represented as

Z̈ = a =
FCφCθ −fZ −mg

m
(4.2)



1034 W.Wang et al.

where Z is the altitude, a is the accelaration of z axis, Cx = cosx, φ, θ are Roll and Pitch
angles, respectively, fZ is the air-resistance, m is themass of theplant and g is the gravitational
acceleration. In general, the Quad-RotorMAVmainlymoves at a low speed, and the lean angle
is about 20 degree. Ignoring the air-resistance, the linearization of Eq. (4.2) can be given

Z̈ = a =
F −mg

m
(4.3)

FromEq. (4.1) and Eq. (4.3), we can get the transfer function from control input to position as
follows

G(s)=
Z(s)

u(s)
=

k

(1+T1s)ms2
(4.4)

To determine the unknown parameters, an information set of input and output data from the
real flight is applied to the system identification tool. In order to verify the accuracy of the
identified model, we analyze the system time historical response, and the result is shown in
Fig. 5. The axis y represents the altitude of theMAV by setting the start point as zero.

Fig. 5. Cross validation result of the altitude model

4.2. Optimal controller design

The altitude model is insensitive to the outer interference. To reduce the control difficulty,
the LQI control method (Wang et al., 2008) is applied. LQI is a preferable design method of
servo controllers based on the optimal control theory.
Rewrite transfer function Eq. (4.4) to the state-space equation

ẋ(t)=Ax+Bu y=Cx (4.5)

To enhance the controller performance in tracing, a newstate variable xr(t) is constructed as the
integral value of the displacement error: xr(t) =

∫

[y(t)−r(t)] dt, and the following expansion
system is established

ẋa(t)=

[

ẋ(t)
ẋr(t)

]

=

[

A 0

−C 0

][

x(t)
xr(t)

]

+

[

B

0

]

u(t)+

[

0

I

]

r(t)=A∗xa(t)+B
∗
u(t)+

[

0

I

]

r(t)

(4.6)

where xa(t) is the state variable of the expanded system, r(t) is the reference input, I is a
identify matrix and

A
∗ =

[

A 0

−C 0

]

B
∗ =

[

B

0

]

(4.7)



Control system design for multi-rotor MAV 1035

The LQI control corresponds to the optimal control problem for the system above in terms of
the quadratic form criterion function

J=
1

2

∞
∫

0

[

x
T
aQxa(t)+u

T(t)Ru(t)
]

dt (4.8)

The optimal control is to determine the control input u(t) tominimize the function. In view of
Eq. (4.8), we have the Hamiltonian

H=
1

2
x
T
a(t)Qxa(t)+

1

2
u
T(t)Ru(t)+λT(t)

[

A
∗(t)xa(t)+B

∗(t)u(t)
]

(4.9)

Then, the adjoint equation and boundary condition can be obtained

λ̇(t)=−
∂H

∂xa
=−A∗

T
(t)λ(t)−Qxa(t) λ(T)=Sxa(T)

∂H

∂u
=Ru(t)+B∗

T
(t)λ(t)=0

(4.10)

with the condition that R is positive, the control input is obtained

u(t)=−R−1B∗
T
(t)λ(t) (4.11)

Assume that the relationship between λ(t) and state variables xa(t) can be expressed as

λ(t)=Pxa(t) (4.12)

in which the matrix P can be obtained from the Riccati equation

PA
∗+A∗

T
P−PB∗R−1B∗

T
P+Q=0 (4.13)

Replacing the λ(t) in Eq. (4.11) by Eq. (4.12), the feedback gain F can be given as follows and
the control input is obtained accordingly

F= [F1,F2] =R
−1
B
∗T
P

u(t)=−F1x(t)−F2xr(t)=−F1x(t)−F2

t
∫

0

e(t) dt
(4.14)

Since the controller is based on the state feedback, it is necessary to observe all the state
conditions. In the altitude model obtained above, there are three state variables. Although the
acceleration and the altitude can be measured by the onboard sensors, there is still one state
variable that cannot be observed. Therefore, a Kalman filter is chosen to estimate the unknown
state variable and it can also filter the noise meanwhile, so that the control quality is improved.

5. Experiments and experiments results

In order to verify the performance of our platform and the technique of attitude and altitude
controller design, real actual flight experiments are carried out. The reference is given to the
MAV by the operator using 72MHz PPM Propo. The flight data are gathered via XBee-pros
which conduct the communication between the aircraft and the ground station.



1036 W.Wang et al.

5.1. Attitude experiments

The effectiveness of the attitude controller is confirmed by keeping the MAV hovering or
moving at a high speed. So the attitude control experiment is divided into two parts: hovering
and tracking experiments.

Results from a hovering test are shown in Fig. 6.We can see from the figure that, when the
reference is set to zero, the airframe can almost keep hovering with deviation of ±3 degrees.
Considering the influences of wind and other interference, the result is acceptable.

Fig. 6. Hovering experimental result

The tracking performance of the controller is shown in Fig. 7. It is clear that the output
angle successfully tracks the reference at both low and high frequencies, indicating that fairly
good tracking results are maintained.

Fig. 7. Attitude experimental result

As to the system robustness which is one of the most interesting features of memory con-
troller, a comparison is conducted between the LQG and memory based control technique for
attitude control with a 300g payload. But once the LQG controller is installed, the flight per-
formance degrades badly, making the operation of the vehicle difficult. Therefore, only the
experimental result for the robustness of memory controller is shown here. Figure 8 indicates



Control system design for multi-rotor MAV 1037

that the controller is almost completely insensitive to model changes and the good robustness
is obtained.

Fig. 8. Robust experimental result

5.2. Altitude experiments

To verify the stability and tracking performance of the designed altitude controller, a step
signal of 5m is given to the Quad-Rotor after hovering stably for some time in a fixed altitude.
Results from the flight test are shown in Fig. 9. The real time experimental result follows the
simulation result well and the output successfully tracks the reference input in about 5s with
a 0.5m overshoot. The controller has also a good stability with the stable control input.

Fig. 9. Altitude experimental result

6. Conclusion

In this research, abrief introduction is given to theQuad-Rotorwhich is selected asour controlled
object, and then the techniques of attitude and altitude controller design are discussed. For
the attitude, memory based controller is proposed due to its insensitivity to the model errors,
parametric uncertainties and other disturbances.Whereas in the altitude control, amodel based
optimal controller is designed. The desired features of the control algorithm are all well proven
in real-time fight experiments.What ismore, by developing thememory controller, the attitude
controller also obtains excellent robustness.

Our future work will focus on indoor obstacle avoidance which is in high demand of the
flight stability. Therefore, on the present basis of the attitude control technique, a relatively
good performance will be achieved.



1038 W.Wang et al.

Acknowledgement

This workwas supported byNationalOutstandingYouth Science Foundation (No. 61125306),Natio-

nal Natural Science Foundation of Major Research Plan (No. 91016004), The National Natural Science

Foundation of China (Grant No. 61004046).

References

1. Achtelik M., Bachrach A., He R., Prentice S., Roy N., 2009, Stereo vision and laser
odometry for autonomous helicopters in gps-denied indoor environments,Proceedings of the SPIE
Conference on Unmanned Systems Technology XI, Orlando, FL

2. Ascending Technologies GmbH, http://www.asctec.de.

3. Bachrach A., He R., Roy N., 2009, Autonomous flight in unstructed and unknown indoor
environments,Proceedings of EMAV

4. HoffmannG.M., HuangH.,Waslander S.L., TomlinC.J., 2011,Precisionflight control for
a multi-vehicle quadrotor helicopter testbed, Control Engineering Practice, 19, 9, 1023-1036

5. http://www.draganfly.com/

6. KendoulF., 2012,Surveyof advances in guidance,navigation,andcontrol ofunmanned rotorcraft
systems, Journal of Field Robotics, 29, 2, 315-378

7. KendoulF., Fantoni I., LozanoR., 2008,Asymptotic stability of hierarchical inner-outer loop-
based flight controllers,Proceedings of the 17th IFAC World Congress, Seoul, Korea, 1741-1746

8. Kendoul F., Lara D., Fantoni I., Lozano R., 2007, Real-time nonlinear embedded control
for an autonomous Quad-Rotor helicopter, AIAA Journal of Guidance, Control, and Dynamics,
30, 4, 1049-1061

9. Martinez V.M., 2007,Modeling of the flight dynamics of a quadrotor helicopter, A MSc Thesis
in Cranfield University

10. Nonami K., Kendoul F., Suzuki S., Wang W., Nakazawa D., 2010, Autonomous Flying
Robots-Unmanned Aerial Vehicles and Micro Aerial Vehicles, Springer

11. Rudin K., Hua M.D., Ducard G., 2011, A robust controller and its application to quadrotor
helicopters,Proceedings of the 18th IFAC World Congress, Italy

12. Song Y.D., 2006, Memory-based control of nonlinear dynamic systems. Part I – Design and
analysis, Industrial Electronics and Applications, 2006 1ST IEEE Conference on

13. Tayebi A., McGilcary S., 2006, Attitude stabilization of a VTOL quadrotor aircraft, IEEE
Transactions on Control Systems Technology, 14, 3, 562-571

14. WangW.,NonamiK.,HirataM., 2008,Model reference slidingmode controlof small helicopter
X.R.B based on vision, International Journal of Advanced Robotic Systems, 5, 3, 235-242

15. Wang W., Song Y.Z., Nonami K., Cheng Y., Zhou Y., Wang F., 2011, Attitude controller
design for a six-rotor typeMAV,Key Engineering Materials, 480-481, 1155-1160

16. Xia M., Zhang Y., Weng L., Ye X., 2012, Fashion retailing forecasting based on extreme
learning machine with adaptive metrics of inputs,Knowledge-Based Systems, 36, 253-259

Manuscript received November 7, 2012; accepted for print June 7, 2013