INTEGRATION OF LIFE SCIENCES AND ENGINEERING -
HIGHLY INFECTIOUS DISEASES IN POPULATION:
CONTROL SYSTEM ANALYSIS AND SYNTHESIS
A THESIS
SUBMITTED TO THE FACULTY OF
THE GRADUATE SCHOOL OF THE UNIVERSITY OF
MINNESOTA
BY
Srijita Bhattacharjee
IN PARTIAL FULFILLMENT OF THE REQUIREMENTS
FOR THE DEGREE OF
MASTER OF SCIENCE
Desineni Subbaram Naidu, Ph.D, Advisor
Peng Fang, Ph.D, Committee Member
Bethany Kubik, Ph.D, Committee Member
July, 2020
Srijita Bhattacharjee, 2020 Copyright c©
All Rights Reserved
Acknowledgement
I would like to express my utmost gratitude to my advisor, Dr. Desineni Subbaram
Naidu, who introduced me to the concept of Convergence which became my sole
motivation to carry out this thesis. Over the course of two years at the University
of Minnesota Duluth (UMD), he has always been a very enthusiastic and inspiring
mentor to me and I truly thank him for his exceptional guidance.
Secondly, I’m expressing my appreciation and gratitude to Dr. Peng Fang and
Dr. Bethany Kubik for serving as committee members of this thesis. Their sincere
advice and corrections have made this work more presentable.
I would like to thank the Electrical Engineering Department of UMD for allowing
me to pursue my MSEE. As a graduate student, I had the opportunity of serving
as a TA in several courses which turned out to be an invaluable experience for me.
It will definitely be very helpful as I proceed with my career ahead.
Lastly, I would like to thank my parents, family and friends for their unceasing
support.
i
Contents
List of Tables iv
List of Figures v
1 Introduction 2
1.1 Literature Survey . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2 Statement of the Problem . . . . . . . . . . . . . . . . . . . . . . . 5
1.3 Overview of the Thesis . . . . . . . . . . . . . . . . . . . . . . . . . 6
2 Nonlinear Modeling of Smallpox in Population 7
2.1 Highly Mobile Population . . . . . . . . . . . . . . . . . . . . . . . 7
2.2 Nonlinear Model of the Disease . . . . . . . . . . . . . . . . . . . . 8
2.3 Simulation of the Nonlinear Model . . . . . . . . . . . . . . . . . . 12
2.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3 Linearization and Control System Characteristics 16
3.1 Linearization of the Nonlinear Model . . . . . . . . . . . . . . . . . 16
3.2 Control System Characteristics . . . . . . . . . . . . . . . . . . . . 18
3.3 Controllability and Observability . . . . . . . . . . . . . . . . . . . 19
ii
3.4 Stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
4 Optimal Control: Regulation Strategy 24
4.1 Optimal Control . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
4.2 Finite-Horizon Linear Quadratic Regulation . . . . . . . . . . . . . 25
4.2.1 Statement of the Problem . . . . . . . . . . . . . . . . . . . 26
4.2.2 Solution for Linear Quadratic Regulation Problem . . . . . . 27
4.3 Simulation Results . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
4.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
5 Conclusion and Future Work 33
5.1 Discussion and Conclusion . . . . . . . . . . . . . . . . . . . . . . . 33
5.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
Bibliography 35
Appendix 38
iii
List of Tables
2.1 New York City Based Parameters . . . . . . . . . . . . . . . . . . . 12
2.2 Initial conditions for nonlinear Smallpox model simulation . . . . . 12
iv
List of Figures
2.1 Simulink Model of the disease . . . . . . . . . . . . . . . . . . . . . 10
2.2 State Response Plots for the Nonlinear Smallpox Model . . . . . . . 14
3.1 Typical Transient Response Characteristics . . . . . . . . . . . . . . 19
3.2 Focus Stability of the Equilibrium State . . . . . . . . . . . . . . . 21
3.3 Knot Stability of the Equilibrium State . . . . . . . . . . . . . . . . 22
4.1 Optimal Control System . . . . . . . . . . . . . . . . . . . . . . . . 28
4.2 Performance of Finite-Horizon Regulation with first set of Q and R 30
4.3 Performance of Finite-Horizon Regulation with second set of Q and R 31
v
Abstract
The research in this master’s thesis focuses on the integration of engineering and
biological sciences. The thesis starts with a literature survey to find out the model
of a highly infectious disease such as smallpox among a certain metropolitan popu-
lation. The chosen nonlinear model is simulated using MATLABR© and SimulinkR©
to obtain input-output relations in the time domain. Next, the nonlinear model
is linearized around an equilibrium point to demonstrate the analysis of the linear
system exhibiting control system characteristics such as controllability, observabil-
ity and stability. Finally, from the control system synthesis (design) point of view,
the theory of linear optimal control has been implemented on the linearized model.
The problem of finite-horizon linear quadratic regulation has been investigated to
provide the necessary optimal feedback to the smallpox model.
1
Chapter 1
Introduction
Convergence of different fields of science in order to come up with new concepts,
technologies as well as software, has been a major area of interest over the past years.
In the history of biomedical revolutions, it is rightfully called the third revolution,
molecular and cellular biology, and genomics being the first two respectively. In
this thesis, a similar concept of converging life science with engineering has been
implemented. The biological model being that of a highly infectious disease like
smallpox analyzed based on the population of a metropolitan city. Convergence
comes into play the moment a biological system is characterized as a mathematical
model. To further explore this integration of life science and engineering, linear
optimal control theory has been implemented on the chosen system.
1.1 Literature Survey
This research was inspired by the theme of Convergence, Integration, and Collabora-
tion Igniting Innovation - Life Sciences (Medicine), Engineering, Arts (Humanities),
2
and Physical Sciences (LEAP) [13, 18, 17, 16]. Convergence has been considered
the third Biomedical Revolution. The first one was Molecular and Cellular Biology
and the second one was Genomics. When two or more fields of science, and in this
case, one being biology and the other being engineering, collaborate with each other,
it develops new perspectives. New concepts of biology, new technologies to study
those concepts, as well as new software to support such a collaboration, are created.
This paves the path for innovation which is a necessity in order to progress. One
of the most fascinating findings of Convergence is the area of precision medicine
which basically refers to diagnostic, and treatments of several diseases. This paved
the path for focusing on specific areas of Biological science and hence was narrowed
down to infectious diseases. Infectious diseases are caused by pathogens and some
of them have the capacity to endanger the life of an individual unless treated prop-
erly. Most of them have the capability of passing on from on person to another.
Unlike a few of the infectious diseases, there are some, which cannot be cured but
only prevented by use of vaccines. Another important factor to keep in mind when
discussing infectious diseases is the probability of mixing among population as that
influences the entire dynamics. Sexually transmitted diseases like AIDS have been
studied under such influences [4]. A highly infectious disease that also is influenced
by such mixing probabilities of people and is air-borne is smallpox which has two
known forms, variola major and minor, and can be transmitted from one individual
to another within an 8 feet radius. Smallpox is characterised by high fever and vom-
iting followed by an outbreak of rash turning into pustules all over the skin. There
always is a possibility that smallpox can be released inspite of its eradication [9].
As stated in [5], smallpox was eradicated in 1980 but two official repositories were
stored, one at the CDC in Atlanta (U.S.) and the second at the State Research
Center of Virology and Biotechnology in Koltsovo, Novosibirsk (Russia). Hence,
the disease can make a comeback in case of some kind of accident or a biological
3
warfare. In an unfortunate scenario of such a bio-terrorism, the vaccination along
with a quarantine campaign is required to stop the outbreak. A daily quarantine
rate of 25% is needed to make that happen [14]. Smallpox is one such disease that
does not have a treatment or cure but can only be prevented by the implementation
of vacciantion. Hence the importance of vaccination is utmost in case of this dis-
ease. The results of [8] indicate that if people voluntarily start to take the vaccine
especially in cases of the first cases, the efficiency of the vaccine in preventing the
outbreak is increased. There also have been studies that indicated to the effects
of antiviral immunity and its duration in an individual after the smallpox vacci-
nation. According to [9] more than 90% of volunteers vaccinated 25-75 years ago
still maintain substantial humoral or cellular immunity (or both) against vaccinia,
the virus used to vaccinate against smallpox. In [11] it has been demonstrated how
population dynamics affect the spread of infectious diseases. But the classification
done there was based on gender while the model selection for this work was done
keeping in mind how the disease could affect a certain population which is highly
mobile increasing the risk of spread. Thus the selection of the model with subway
and nonsubway users was chosen. The basic theory used to model the disease in [6]
included three stages of infection i.e. susceptible, infected, and recovered. The stage
that comes before being infected is exposure to the disease and hence the model
used in this work has considered the mobility of people as an important factor as
that influences spread of a disease to a great deal. Hence the entire population of
a metropolitan city like New York has been divided into subway and non-subway
users. They again have been divided into four categories of susceptible, exposed,
infectious, and recovered individuals [3]. A similar model was also used in [12] where
they have established immunosuppresion to be stronger than vaccination but in this
work the focus is mainly on the vaccination rate. The nonlinear modeling has been
done and the equilibrium point has been carefully chosen in accordance with the
4
behavior of the state responses [10]. The parameters for this work have been taken
from [5]. The concept that a system may or may not be controllable was intro-
duced by Kalman and has been explored in this work. The necessary conditions for
positive controllability as in [22] has been tested and discussed. Similarly, another
control system characteristic that has been explored in this work is the observability
of the system. As stated in [20], the rank of the observability matrix is derived to
discuss it later. The Lyapunov equation can be solved using the system matrix A
to derive P as shown in [23].
ATP + PA+ I = 0 (1.1.1)
where I is identity matrix. By judging if P is positive definite or not the stability
of the system can be judged. In this work the stability of the equilibrium point
has been explored by referring to the deductions mentioned in [19]. Hence the
eigenvalues of the system matrix have been observed to come to a conclusion. In
[7] it has been demonstrated how sensitivity analysis might help in understanding
and control a smallpox outbreak. The final section of this work is to investigate
the finite-horizon linear quadratic regulation in order to derive a control method for
bringing down the susceptibility. As established in [2] and [21], the choice or tuning
of the weighted matrices Q and R determines the final regulation result.
1.2 Statement of the Problem
Biological systems can be studied using engineering tools and techniques which can
give rise to new concepts, technologies, and software. This initiative, also known
as convergence, the third revolution in the history of biomedical revolutions, can
provide advanced control strategies for such biological systems by studying and
5
analyzing their behavior. This work focuses on implementation of linear optimal
control theory. It involves linearization of such a nonlinear system around a carefully
chosen equilibrium point so that the control system characteristics of it can be
studied. Further, the finite-horizon linear quadratic regulation problem has been
addressed.
1.3 Overview of the Thesis
This dissertation is composed of five chapters covering the following topics:
Chapter 2 provides a detail discussion on the modeling of the infectious disease,
smallpox in population. In particular, discussions on how the population of New
York City is divided into certain groups of people in order to identify them in terms
of the disease is provided here. Additionally, the nonlinear model is simulated and
the responses are analyzed in this section. Chapter 3 presents the discussion on
linearization of the nonlinear model of smallpox around the equilibrium point. It
also deals with the control system characteristics like controllability, observability,
and stability of the linearized system. Chapter 4 starts with a discussion on optimal
control and then discusses the finite-horizon linear quadratic regulation problem.
The solution is examined using the matrix differential Riccati equation and the
simulation results are studied. Finally, Chapter 5 presents the conclusion of the
work with a note on future improvement scopes.
6
Chapter 2
Nonlinear Modeling of Smallpox
in Population
2.1 Highly Mobile Population
The infectious disease, smallpox, was modeled based on a certain population of a
metropolitan city. The city that was modeled was divided into N neighborhoods.
The population within each of those neighborhood was further divided into two
group of individuals- the subway users (SU) and the nonsubway users (NSU). The
SU individuals were assumed to have been in contact with both SU individuals all
over the city and NSU individuals but only within their home neighborhood. Si-
multaneously, the NSU individuals were assumed to have been in contact with only
individuals living in the same neighborhood, since they did not access the subway.
This type of modeling was done based on the assumption that the individuals here
did not change their mass-transportation status. In case of a smallpox attack, not
7
all individuals would maintain the same epidemiological status. Hence, they are fur-
ther divided into specific epidemiological classes of susceptible, exposed, infectious,
and recovered denoted by Wi, Xi, Yi, Zi for the SU individuals and Si, Ei, Ii, Ri for
the NSU individuals respectively. The subscript of i denotes the ith neighborhood.
In this work, only the NSU individuals have been focused upon. Thus, the
model presented in the following section solely focuses on the individuals off the
subway and for one particular neighborhood, hence i = 1. The data considered here
is that of New York City in 2000.
2.2 Nonlinear Model of the Disease
A fourth order, nonlinear, state-space model of smallpox is presented. The dynamic
equations are
dS
dt
= Re−B − (µ+ ql1)S,
dE
dt
= B − (µ+ φ+ ql2)E,
dI
dt
= φE − (µ+ φ+ d)I,
dR
dt
= αI − µR + ql1S + ql2E,
(2.2.1)
where, S,E, I, and R denote the susceptible, exposed, infectious, and recovered
individuals. Re is considered to be the recruitment rate of the individuals assumed
to be 0.7 which includes recruitment both by birth and immigration. The NSU
individuals are infected at a rate of B. The per capita natural mortality rate is µ
and that due to smallpox specifically is d. l1 and l2 are the vaccination efficacy rates
of susceptible and exposed individuals respectively. The per capita progression rate
from latent to the infectious is denoted as φ. The rate at which infected individuals
recover is α. The number of susceptible individuals is calculated by considering the
8
people responding to the vaccine as well people dying in the neighborhood. These
two terms along with the infection rate B is deducted from the recruitment rate Re
so the remaining is all the susceptible individuals in the first equation. The rate φ
with which the exposed individuals progress from latent to the infectious state is
subtracted from the second equation and added to the third, since once infected they
are no more just exposed to the disease but also infected by it. All the individuals
who have come in direct or indirect contact with the virus are considered exposed
to the disease. They do not include the ones dead, already infected, and responding
to the vaccine and hence the term (µ+φ+ ql2) is subtracted from the infection rate
B(t) in the second equation. Individuals dying at the rate of µ, dying specifically
due to smallpox at a rate of d, and those still progressing towards more infection at
the rate of φ are eliminated from all the exposed people who have progressed to be
infected at the rate of φ in the third equation. The recovered individuals depend
on all the others. Recovered are the ones who were susceptible to the disease but
responded to the vaccine (ql1), the ones exposed to the disease but still responded
to the vaccine (ql2) and the infected individuals recovering at the rate of α. The
ones eliminated from the fourth equation are the individuals who pass away at the
rate of µ. The input here is the vaccination rate q which has been assumed to be a
step function. Keeping it at a moderate level aids in analysing the original system.
The vaccination rate cannot be assumed to be 1 here because practically that is
impossible to achieve. Hence, from the modeling of the disease it is quite clear
that infection from the disease can happen to an individual only if exposed to the
disease. Exposure could be due to direct contact with another infected individual
or indirect contact with their used items. Recovery does not only indicate recovery
from the disease itself at a natural rate once infected but also by vaccinating those
susceptible or exposed to the disease.
The nonlinear model is built using SimulinkR© as in figure 2.1 shown below. The
9
block named as subsystem expands further to hold the functions forming the four
states.
Figure 2.1: Simulink Model of the disease
The infection rate B is represented as-
B = βaS
[
PaI
Tτ +Q
+
PbY τ
Tτ +Q
]
, (2.2.2)
where β is the transmission rate per contact, a is the average number of contacts
of NSU per day, Y is the number of infected SU individuals assumed to be 1
million individuals, Pa is the mixing probability of SU individuals, Pb is the mixing
probability of SU and NSU individuals which means
Pa + Pb = 1. (2.2.3)
Hence (2.2.2) comes down to
B =
βaS
Tτ +Q
[(1− Pb)I + (PbY τ)] . (2.2.4)
10
Here, T is the total number of SU population taken to be 0.2 million per day in
this study, τ is the proportion of time spent by SU individuals off the subway, Q is
the total NSU population which is approximately 8 million here, Y is the infected
SU population and Pb is expressed as
Pb =
bτT
aQ+ bτT
, (2.2.5)
where a and b are the average number of contacts of NSU and SU per day given
as 5 and 10 respectively and τ is expressed as
τ =
σ
σ + ρ
, (2.2.6)
where ρ and σ denote the per capita rates of SU individuals at which they get
on and off the subway. We have assumed an approximate value for τ to be 0.6 as
that is a standard proportion of time of SU individuals spent off the subway. The
infection rate B is generated as 0.152 using the parameters given. Though it is only
the NSU individuals being discussed here but to calculate the infection rate, the
total infectious SU individuals had to be assumed to a standard value because they
bring in the infection from the all the travelling and using of subway. That has a
huge part to play in infecting the NSU individuals.
The numerical values of the constant parameter of this model are listed in Table
2.1
11
Table 2.1: New York City Based Parameters
Parameter Description Value
beta transmission rate per contact 0.5
a Per capita average no. of contacts of NSU per unit of time 5
b Per capita average no. of contacts of SU per unit of time 10
d Per capita mortality rate due to smallpox 0.0116
α recovery rate 0.086
µ Per capita natural mortality rate 0.033
l1 Vaccination efficacy rate for susceptible population(S) 0.97
l2 Vaccination efficacy rate for exposed population(E) 0.3
φ Per capita infection rate from latent to infectious 0.086
Q Total NSU population 8 ∗ 106
T Total SU population 0.2 ∗ 106 per day
2.3 Simulation of the Nonlinear Model
To gain insights into its dynamic behavior, the nonlinear Smallpox model (2.2.1)
is simulated with the numerical parameters in Table 2.1. The initial conditions
for the state variables are set as in Table 2.2. The simulations were performed in
MATLABR© and SimulinkR© software.
Table 2.2: Initial conditions for nonlinear Smallpox model simulation
Parameter Description Value
S(0) Susceptible Individuals 10
E(0) Exposed Individuals 8
I(0) Infectious Individuals 0
R(0) Recovered Individuals 0
The response of the states to the initial conditions is provided in Figure 2.2.
Initially the number of infectious and recovered individuals are considered none.
Number of susceptible individuals is assumed to be 10 and exposed individuals to
be 8 in each neighborhood. Susceptibility to the disease can be due to age, poor
immunity, no vaccination, pregnancy and other medical conditions which restricts
them from taking the vaccine if needed. 8 exposed individuals is assumed based on
12
the fact that they might have been directly or indirectly in contact with the infected
SU individuals (Y ). The equilibrium points are calculated by setting the equations
in (2.2.1) to zero and solving them as in (2.3.1),
Re−B(t)− (µ+ ql1)S∗ = 0,
B(t)− (µ+ φ+ ql2)E∗ = 0,
φE∗ − (µ+ φ+ d)I∗ = 0,
αI∗ − µR∗ + ql1S∗ + ql2E∗ = 0,
(2.3.1)
The values are rounded up because they represent number of people.
(S∗, E∗, I∗, R∗) = (2, 1, 1, 19) (2.3.2)
The response plots demonstrate that the recovery of people depend more on the
susceptible and exposed individuals and their response to the vaccine than the ones
already infected. Around 2 among the 8 exposed have a chance of getting infected
while most of the remaining population respond to the vaccine. But once infected,
it takes approximately a duration of four weeks for them to recover from it if at
all they do. The number of recovered individual surges to 11 in 4 days while the
number of exposed people drops to 3 from 8 in 4 days as that is the standard time
period within which if the vaccine is provided then the infection can still be stopped.
However, the drop in number of susceptible population within 3 days is from 10 to
3 which is higher than the exposed population because of the higher vaccination
efficacy rate (l1) of the susceptible individuals. They never go to a value of zero
which is a cause of concern since that indicates that susceptibility to the disease
can always be present in the population due to several reasons like pregnancy, old
age, no vaccination etc. The recovered individuals also never go to a straight zero
13
(a) S(t) (b) E(t)
(c) E(t) (d) R(t)
Figure 2.2: State Response Plots for the Nonlinear Smallpox Model
because of the same reason. All of them reach the equilibrium point at around a
month’s time which is 30 days. This makes sense because if and when infected by
the virus, that is the exact amount of time needed to recover from the disease if
that is a possibility.
2.4 Conclusion
In this chapter the nonlinear model of the smallpox disease has been presented along
with a discussion on how that has been done in a certain population of New York
14
City. The parameters have been explained in order to explain the modeling. The
nonlinear model has been built and simulated using MatlabR© and SimulinkR© and
the responses have been analyzed within a particular time period.
15
Chapter 3
Linearization and Control System
Characteristics
3.1 Linearization of the Nonlinear Model
The nonlinear smallpox model was linearized around the equilibrium point and the
eigenvalue was evaluated. The linearization was carried out in MATLABR© and
SimulinkR©. The values from Table 2.1 were used for the simulation.
In order to study the behavior of a nonlinear dynamical system, we usually
linearize it around an operating point. In this case, the time instant at which it
is linearized is 30 days. The operating point at that time instant is (1.6, 0.1, 0.2,
18). 30 days is crucial for a disease like smallpox because if an individual exposed
to the disease gets infected inspite of the vaccination, then with an approximately
10% recovery rate the infectious individuals have a chance of completely recovering
after four weeks. At the end of four weeks when all the scabs from the rashes have
16
fallen off, the person is not contagious anymore. Hence 30 days is the time duration
required for the system to become stable.
Suppose, the nonlinear system is,
x˙(t) = f(x) + g(x)u(t),
y(t) = h(x),
(3.1.1)
where, x is nth state vector, u(t) is rth control vector and y(t) is mth output
vector. The system is transformed to a linear structure as
x˙(t) = Ax(t) + Bu(t),
y(t) = Cx(t),
(3.1.2)
where, f(x) = Ax(t), B= g(x), and h(x) = Cx(t). Here, A is n×n state matrix
and B is n× r control matrix.
After linearizing around the operating point of (1.6, 0.1, 0.2, 18) the system that
is obtained is given as,
A =

−0.4369 0 0 0
0.0159 −0.2390 0 0
0 0.0860 −0.1306 0
0.3880 0.1200 0.0860 −0.0330
 ;B =

−1.554
−0.032
0
1.586

C =

1 0 0 0
0 1 0 0
0 0 1 0
0 0 0 1
 ;D =

0
0
0
0

(3.1.3)
The eigenvalue obtained for the linearized system matrix is
17
E =

−0.033
−0.2390
−0.1306
−0.4369
 (3.1.4)
The importance of linearization of the given nonlinear model is to be able to
study the system characteristics as in the following sections. Further, it paves
the way for solving the linear quadratic regulation problem. The eigenvalues are
important because they help us denote the stability of the equilibrium point.
3.2 Control System Characteristics
The desired performance characteristics of control systems are practically described
in terms of time domain quantities usually. The response is generally observed by
applying a unit step input because if that is known, then mathematically it is pos-
sible to compute for any other input. Some of the transient response characteristics
most commonly used to analyze the system are as follows,
• Delay Time (Td) - The time taken by the response to reach half of its final
value initially.
• Rise Time (Tr) - The time taken by a system response to reach a specified
high value from a specified low value.
• Peak Time (Tp) - The time taken for the response to reach the first peak of
the overshoot.
• Maximum Overshoot (Mp) - It is the maximum peak value of the response
curve calculated from unity.
18
• Settling Time (Ts) - It is defined as the time required for the response to reach
and stay within a range about the final value of size specified by the absolute
percentage of the final value.
Figure 3.1: Typical Transient Response Characteristics
img src: https://electronicsguide4u.com/control-system-transient-response-specification/
The Figure 3.1 depicts all the above mentioned characteristics.
3.3 Controllability and Observability
Controllability can be defined as the ability of an external input to change the initial
internal state of a system to any desired state in a finite time interval. There are
more than one condition that can determine if a system is controllable. In this
study, one of them has been utilized which is to find the rank of the controllability
19
matrix given by
Qc =
[
B AB A2B . . . An−1B
]
(3.3.1)
The rank denotes the number of linearly independent rows and columns of the
controllability matrix Qc. If it equals a full rank n meaning if the rank is equal to
the order of the system, then we conclude that the system is controllable.
The linearized smallpox system (3.1.3) is subjected to a controllability test Qc
is derived as
Qc =

−1.5540 0.6789 −0.2966 0.1296
−0.0320 −0.0171 0.0149 −0.0083
0 −0.0028 −0.0011 0.0014
1.5860 −0.6591 0.2829 −0.1227
 (3.3.2)
which indicates 4 linearly independent rows and columns. Hence the rank of the
controllability matrix is 4 which is a full rank which means the system is controllable.
Thus, the control input, which is the vaccination rate (q) in this case, has the ability
to change the initial states to any other state which is what happens practically.
Smallpox is one of those diseases which do not have a medicated treatment but is
solely dependent on the vaccine and thus the system is rightfully controllable by its
input.
A given system is said to be observable if it is possible to determine the initial
states of the system by observing the outputs in finite time interval. There are more
than one condition that can determine if a system is observable. In this study, one
of them has been utilized which is to find the rank of the observability matrix given
by
Qo =
[
C CA CA2 . . . CAn−1
]T
(3.3.3)
The rank denotes the number of linearly independent rows and columns of the
observability matrix Qo. If it equals a full rank n meaning if the rank is equal to
20
the order of the system, then we conclude that the system is observable.
The linearized smallpox system (3.1.3) is subjected to an observability test and
Qo is tested to have a rank of 4 which is a full rank which means the system is This
indicates that it is very well possible to gauge the initial states of the system by
studying the output response of each of them at different time intervals.
3.4 Stability
Stability of a dynamical system denotes that the trajectories do not deviate much
under small perturbation. Hence, perturbing the initial state in a particular di-
rection results in the trajectory asymptotically approaching it and in case of other
directions to the trajectory getting away from it. If the eigenvalues have an imag-
inary part then the focus represents the stability as shown in figure 3.2 where the
numbers 1, 2, 3, and 4 represent stable, asymptotically stable, globally asymptoti-
cally stable, and unstable equilibrium state respectively.
Figure 3.2: Focus Stability of the Equilibrium State
However, when the eigenvalues have no imaginary part, then the stability is
represented by a knot as shown in Figure 3.3
Asymptotic stability refers to whether a nearby orbit of a trajectory converges
21
Figure 3.3: Knot Stability of the Equilibrium State
img src: https://web.ma.utexas.edu/users/davis/375/popecol/lec9/equilib.html
with the given orbit. The qualitative behavior of the said orbit under perturbations
can be studied by linearizing the system around an equilibrium point which is also
the centre for the orbit. The eigenvalues of the system matrix A determines the
stability of the system by characterizing the behavior of the nearby operating points.
If all eigenvalues are negative real numbers then the system will be driven back to
its steady-state when disturbed. The eigenvalues for the system matrix
A =

−0.4369 0 0 0
0.0159 −0.2390 0 0
0 0.0860 −0.1306 0
0.3880 0.1200 0.0860 −0.0330
 is E =

−0.033
−0.2390
−0.1306
−0.4369

All the eigenvalues are negative real numbers indicating that the system (3.1.3)
is asymptotically stable. This means that the equilibrium point around which the
model is linearized is a stable one. Hence the responses should approach zero as
22
time approaches infinity. This is true in this case where the number of exposed and
infectious individuals reduces to zero as time increases making it an asymptotically
stable system.
3.5 Conclusion
In this chapter, the nonlinear smallpox model (2.2.1) has been linearized around an
equilibrium point which is chosen based on the recovery of the individuals already
infected with the disease. The eigenvalues of the system matrix is derived. In the
next section, three of the control system characteristics, controllability, observabil-
ity, and stability, have been discussed.
23
Chapter 4
Optimal Control: Regulation
Strategy
4.1 Optimal Control
Optimal control aims to optimize a plant, in this case a process or system, by finding
a control input. To obtain a control input, the optimal control minimizes or max-
imizes a process by satisfying a specific performance criterion. The system should
have the following three characteristics in order to be considered as an optimal
control problem:
1. A mathematical model of the system that needs to be controlled. The model
is generally expressed as a dynamic system with state variables.
2. Description of a performance index such as reaching a target in a minimal
amount of time.
24
3. Description of boundary conditions or physical constraints needed to exceed
to reach the target.
An optimal control problem is generally derived by Pontryagin’s minimum principal
or by Hamilton-Jacobi-Bellman (HJB) equation. Further, calculus of variation plays
an important role in obtaining a solution for optimal control problems. Examples
of optimal control problems are
1. Reaching moon with minimum fuel expenditure ( Fuel- Optimal Control Sys-
tem). Mathematical model for this system would be a dynamic system of a
spacecraft with associated state variables such velocity, thrust of a rocket.
2. Reaching from point A to point B in minimal time, etc.
3. Some of the areas that can be optimally controlled are biomedical systems,
aerospace systems, automotive systems etc.
Here, linear control theory has been applied to solve the linear quadratic reg-
ulator problem. By implementing linear control theory, state regulation has been
achieved.
4.2 Finite-Horizon Linear Quadratic Regulation
Linear quadratic regulation (LQR) is the procedure to regulate and keep the state
near zero at any given interval of time. This can be necessary in various scenarios
for different practical systems. Those which respond to sudden disturbances due
sudden change in any system parameter could be subjected to such a regulation.
Hence, the intention is to bring the state from a non-zero value to a zero value in
the finite amount of time. In this work the intention is to bring one of the states to
zero, hence it is called a state regulator problem.
25
In the following section, all the important elements required to achieve linear
quadratic regulation are described.
4.2.1 Statement of the Problem
The system (2.2.1) has been linearized to obtain a structure which is given as,
x˙(t) = Ax(t) +Bu(t), (4.2.1)
y(t) = Cx(t), (4.2.2)
Here, A is n × n state matrix and B is n × r control matrix. And, x(t) is the nth
order state vector, u(t) is the rth order control vector, and y(t) is the mth order
output vector. Let z(t) be the mth order desired output.
and The performance index is chosen to be
J =
1
2
x′(tf )F (tf )x(tf ) +
1
2
∫ tf
t0
[x′(t)Q(t)x(t) + u′(t)R(t)u(t)]dt (4.2.3)
with tf specified and x(tf ) not specified which indicates that it is a free-final
state system. Here u(t) is not constrained. Also it is assumed that F (tf ) and Q(t)
are n× n symmetric, positive semi-definite matrices, and R(t) is r × r symmetric,
positive definite matrix.
26
4.2.2 Solution for Linear Quadratic Regulation Problem
The optimal control u∗(t) is given by,
u∗(t) = −R−1(t)B′(t)P (t)x∗(t)
= −K(t)x∗(t)
(4.2.4)
where
K(t) = R−1(t)B′(t)P (t) (4.2.5)
is called Kalman gain and P (t), the n × n symmetric, positive definite matrix
(∀ t ∈ [t0, tf ]), is the solution of the matrix differential Riccati equation (DRE)
P (˙t) = −P (t)A(t)− A′(t)P (t)−Q(t) + P (t)B(t)R−1(t)B′(t)P (t) (4.2.6)
with final condition
P (t = tf ) = F (tf ) (4.2.7)
The optimal state is the solution of the linear state equation
x˙∗(t) = [A(t)−B(t)R−1(t)B′(t)P (t)]x∗(t) (4.2.8)
and the optimal cost J∗ where
J∗ =
1
2
x∗
′
(t)P (t)x∗(t) (4.2.9)
27
The optimal control u∗(t), given by (4.2.4), is linear in the optimal state x∗(t).
Figure 4.1: Optimal Control System
Implementation of the optimal control system is shown in the figure 4.1[15].
28
4.3 Simulation Results
This section focuses on the simulation results. For the following linear system (4.3.1)
A =

−0.4369 0 0 0
0.0159 −0.2390 0 0
0 0.0860 −0.1306 0
0.3880 0.1200 0.0860 −0.0330
 ;B =

−1.554
−0.032
0
1.586

C =

1 0 0 0
0 1 0 0
0 0 1 0
0 0 0 1
 ;D =

0
0
0
0

(4.3.1)
Consider the first set of weighted matrices as
Q =

1000 50 50 50
50 100 50 50
50 50 50 50
50 50 50 50

and
R =
[
30800
]
The value of these weighted matrices are extremely essential in such regulator prob-
lems. Without a proper set of weighted matrix, usually it is difficult to make the
state reduce to a zero value. Therefore, value of Q and R were chosen by trial and
error process to achieve the best regulation case.
The aim here was to achieve a zero tracking for one of the states of the system.
29
The state chosen for regulation was the number of susceptible individuals. We had
observed in the original system that susceptibility was existent to a constant of
2 people even after 60 days of time period. The control of the linear quadratic
regulator is tuned using the weighted matrices Q and R so that the susceptible
individuals track a zero value in minimum number of days which in this case is
about a week. This indicates to the possibility of controlling the susceptibility in
the neighborhood by optimizing the vaccination. It may also indicate to the fact
that susceptibility to the disease is majorly due to lack of vaccination which in turn
may suggest a mass vaccination as compulsory before anybody is actually exposed.
The result is shown in the Figure 4.2.
(a) Regulated susceptible individuals (b) Optimized u
Figure 4.2: Performance of Finite-Horizon Regulation with first set of Q and R
Similarly, the weighted matrices are again tuned to obtain a second set as follows,
Q =

200 0 0 50
0 50 50 50
0 50 50 50
50 50 50 50

30
and
R =
[
5000
]
These set of Q and R gives the result as shown in figure 4.3. The Susceptible
Individuals are regulated to track a zero value but due to change in the tuning
parameters Q an R, it reaches a zero value in about 14 days which is less efficient
than the regulation result in 4.2. Hence, this highlights the importance of efficiently
tuning the weighted matrices as that directly affects the optimization of the control
of the linear quadratic regulator. As the control is modified, so is the state. Thus,
though both the set of Q and R successfully regulate the state but tuning them
efficiently is important as that directly impacts the time period within which the
state is regulated to zero. Hence, comparing the two results, it is quite evident that
the first set of Q and R giving the result found in Figure 4.2 is more efficient.
(a) Regulated susceptible individuals (b) Optimized u
Figure 4.3: Performance of Finite-Horizon Regulation with second set of Q and R
This method of optimizing the control input in bringing down the susceptibility
is required because the adverse effects of the smallpox vaccine does not allow a
mass vaccination. Also, it has already been proved in [1] that mass vaccination
is not possible. Hence, finding the optimized vaccination rate to bring down the
susceptibility to zero has been aimed for here.
31
4.4 Conclusion
This chapter begins with a procedural description of the theory of optimal control
which lays grounds for the finite-horizon linear quadratic regulation problem. The
state of susceptibility is regulated by tuning the weighted matrices to two different
set of values and the results are discussed.
32
Chapter 5
Conclusion and Future Work
5.1 Discussion and Conclusion
The implementation of linear control theory on a nonlinear infectious disease model
was carried out in several steps. The work began with analyzing the nonlinear
model. The next step was to linearize it around an equilibrium point which was
chosen after studying the response of human body to the disease. Further, the
control system characteristics like controllability, observability, and stability of the
system were discussed. Finally, linear optimal control strategy was applied to opti-
mize the state of susceptibility in order to regulate it to a zero value.
A paper titled “Integration of Life Sciences and Engineering- Highly Infectious
Diseases in Population”, which summarizes the initial analysis done in this work,
authored by Srijita Bhattacharjee, and Dr. Desineni Subbaram Naidu, was submit-
ted and accepted for The 24th World Multi-Conference on Systemics, Cybernetics
and Informatics: WMSCI 2020, to be held on September 13 - 16, 2020 virtually.
33
5.2 Future Work
Considering the fact that this work was carried out on a part of the entire model,
the following goals can be accomplished further.
1. Finding an advanced finite-horizon nonlinear tracking strategy for the states.
Therefore, with SDRE technique, a SOC dependent tracking strategy devel-
opment can be implemented to study the results.
2. Implementation of advanced control strategies to lower the possibility of in-
fection after the period of exposure or in other words trying to control the
vaccination efficacy rate of the exposed individuals to reduce infection.
34
Bibliography
[1] GK Aldis and MG Roberts. An integral equation model for the control of a
smallpox outbreak. Mathematical biosciences, 195(1):1–22, 2005.
[2] R. Ali, M. Liaqat, F. M. Malik, and M. B. Malik. Optimal output quadratic
regulation. In 2016 6th International Conference on Intelligent and Advanced
Systems (ICIAS), pages 1–4, 2016.
[3] Fred Brauer and Carlos Castillo-Chavez. Mathematical Models for Communi-
cable Diseases. Society for Industrial and Applied Mathematics, Philadelphia,
PA, USA, 2012.
[4] Stavros Busenberg and Carlos Castillo-Chavez. A general solution of the prob-
lem of mixing of subpopulations and its application to risk-and age-structured
epidemic models for the spread of aids. Mathematical Medicine and Biology:
A Journal of the IMA, 8(1):1–29, 1991.
[5] Carlos Castillo-Chavez, Baojun Song, and Juan Zhang. 8. An Epidemic Model
with Virtual Mass Transportation: The Case of Smallpox in a Large City, pages
173–197. SIAM, Philadelphia, PA, 2003.
35
[6] Neil M. Ferguson, Matt J. Keeling, W. John Edmunds, Raymond Gani,
Bryan T. Grenfell, Roy M. Anderson, and Steve Leach. Planning for smallpox
outbreaks. Nature, 425, 2003.
[7] William H Foege, J Donald Millar, and J Michael Lane. Selective epidemiologic
control in smallpox eradication. American journal of epidemiology, 94(4):311–
315, 1971.
[8] M Elizabeth Halloran, Ira M Longini, Azhar Nizam, and Yang Yang. Contain-
ing bioterrorist smallpox. Science, 298(5597):1428–1432, 2002.
[9] Erika Hammarlund, Matthew W Lewis, Scott G Hansen, Lisa I Strelow, Jay A
Nelson, Gary J Sexton, Jon M Hanifin, and Mark K Slifka. Duration of antiviral
immunity after smallpox vaccination. Nature medicine, 9(9):1131, 2003.
[10] Shaleena Jaison and Desineni Naidu. Integration of life sciences and engineering
- optimal control of hiv using time scales. pages 1–3, 03 2019.
[11] K. Koide, H. Matsuura, N. Noda, and M. Nakano. Simulations on the popu-
lation and disease in japan. In First International Conference on Innovative
Computing, Information and Control - Volume I (ICICIC’06), volume 3, pages
550–552, 2006.
[12] C Raina MacIntyre, Valentina Costantino, Xin Chen, Eva Segelov, Abrar Ah-
mad Chughtai, Anthony Kelleher, Mohana Kunasekaran, and John Michael
Lane. Influence of population immunosuppression and past vaccination on
smallpox reemergence. Emerging infectious diseases, 24(4):646, 2018.
[13] Massachusetts Institute of Technology (MIT), Cambridge, MA, USA. The
Third Revolution: The Convergence of Life Sciences, Physical Sciences and
Engineering, 2011. The First Revolution: Molecular and Cellular Biology and
The Second Revolution: Genomics.
36
[14] Martin I Meltzer, Inger Damon, James W LeDuc, and J Donald Millar. Mod-
eling potential responses to smallpox as a bioterrorist weapon. Emerging in-
fectious diseases, 7(6):959, 2001.
[15] D Subbaram Naidu. Optimal control systems. CRC press, 2002.
[16] The National Academies of Science, Engineering and Medicine, Washington,
DC, USA. Collaborations of Consequence: NAKFI s 15 Years Igniting Inno-
vation at the Intersections of Disciplines, 2018. The National Academies Keck
Futures Initiative (NAKFI.
[17] The National Academies of Science, Engineering and Medicine, Washington,
DC, USA. The Integration of the Humanities and Arts with Sciences, Engi-
neering, and Medicine in Higher Education: Branches from the Same Tree,
2018. Contributed by David Skorton and Ashley Bear, Editors; Committee on
Integrating Higher Education in the Arts, Humanities, Sciences, Engineering,
and Medicine; Board on Higher Education and Workforce; Policy and Global
Affairs; National Academies of Sciences, Engineering, and Medicine.
[18] National Research Council of The National Academies, Washington, DC, USA.
Convergence: Facilitating Transdisciplinary Integration of Life Sciences, Phys-
ical Sciences, Engineering, and Beyond, 2014. Contributed by Committee on
Key Challenge Areas for Convergence and Health; Board on Life Sciences;
Division on Earth and Life Studies; National Research Council.
[19] Marc R Roussel. Stability analysis for odes. Nonlinear Dynamics, lecture notes,
University Hall, Canada, 2005.
[20] K. Shen, A. V. Proletarsky, and K. A. Neusypin. Quantitative analysis of
observability in linear time-varying systems. In 2016 35th Chinese Control
Conference (CCC), pages 44–49, 2016.
37
[21] S. Ullah and M. Liquat. Optimal output regulation on sample-data systems.
In 2015 International Conference on Control, Electronics, Renewable Energy
and Communications (ICCEREC), pages 79–84, 2015.
[22] H. Yoshida and T. Tanaka. Positive controllability test for continuous-time
linear systems. IEEE Transactions on Automatic Control, 52(9):1685–1689,
2007.
[23] X. Zeng, Z. Li, Y. Wu, W. Gao, J. Zhang, M. Ren, and B. Zhang. Dynamic
stability analysis based on state-space model and lyapunov’s stability criterion
for sic-mos and si-igbt switching. In 2018 IEEE 30th International Symposium
on Power Semiconductor Devices and ICs (ISPSD), pages 268–271, 2018.
38
Appendix
Appendix presents MATLAB codes that were used in this thesis. For the electronic
version of these codes please contact Dr. D Subbaram Naidu (dsnaidu@d.umn.edu).
39
Code for the parameters
1 clear a l l
2 close a l l
3 clc
4 % ∗∗∗∗∗∗∗∗∗∗ Thesis Model Parameters ∗∗∗∗∗∗∗∗∗∗∗∗ %
5 Re = 0 . 7 ;
6 alpha = 0 . 0 8 6 ;
7 mu = 0 . 0 3 3 ;
8 d = 0 . 0 1 1 6 ;
9 l 1 = 0 . 9 7 ;
10 l 2 = 0 . 3 ;
11 phi = 0 . 0 8 6 ;
12 beta = 0 . 5 ;
13 a = 10 ;
14 b = 30 ;
15 tau = 0 . 6 ;
16 Y= 1∗10ˆ6;
17 T = 0.2∗10ˆ6 ;%per day
18 Q = 8∗10ˆ6;%per day
19 %S = 3∗10ˆ6; Y = 0;
20 P = [ b∗ tau∗T ] / [ [ a∗Q]+[b∗ tau∗T ] ] ;
21 Z = [T∗ tau ]+Q;
22 % % ∗∗∗∗∗∗∗∗∗∗∗∗∗ I n i t i a lC ond i t i o n s ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗ %
23 i c x 1 = 10 ;
24 i c x 2 = 8 ;
25 i c x 3 = 0 ;
26 i c x 4 = 0 ;
40
27 %% Eqm. Point
28 %S = (Re−B) /(mu+(q∗ l 1 ) ) ;
29 %E = B/(mu+phi+(q∗ l 2 ) ) ;
30 %I = phi ∗E/(mu+alpha+d) ;
31 %R = (( a lpha ∗ In )+(q∗ l 1 ∗Sus )+(q∗ l 2 ∗Ex) )/mu;
32 %Run t h i s f o l l owed by the Simul ink Model
41
Code for linearization
1 %% Fi r s t run the s imu l ink model f o l l owed by t h i s f i l e
2 clc ;
3 %I/O ob j e c t and Op. po in t o b j e c t
4 d i sp l ay ( ’ ∗∗∗∗∗∗∗∗∗∗ I /O ob j e c t c r ea ted ∗∗∗∗∗∗∗∗∗∗∗∗∗ ’ )
5 i o=g e t l i n i o ( ’Newmodel ’ ) ;
6 d i sp l ay ( ’ ∗∗∗∗∗∗∗∗∗Operating Point Object ∗∗∗∗∗∗∗∗∗∗∗∗∗∗ ’ )
7 op po int=f indop ( ’Newmodel ’ , [ 4 7 14 30 6 0 ] ) ;
8 %More than one opera t ing po in t . .
9 d i sp l ay ( ’ ∗∗∗∗∗∗∗∗1 s t op∗∗∗∗∗∗∗ ’ )
10 op po int (1 )
11 d i sp l ay ( ’ ∗∗∗∗∗∗∗∗2nd op∗∗∗∗∗∗∗ ’ )
12 op po int (2 )
13 d i sp l ay ( ’ ∗∗∗∗∗∗∗∗3 rd op∗∗∗∗∗∗∗ ’ )
14 op po int (3 )
15 d i sp l ay ( ’ ∗∗∗∗∗∗∗∗4 th op∗∗∗∗∗∗∗ ’ )
16 op po int (4 )
17 d i sp l ay ( ’ ∗∗∗∗∗∗∗∗5 th op∗∗∗∗∗∗∗ ’ )
18 op po int (5 )
19 %d i s p l a y ( ’∗∗∗∗∗∗∗∗∗∗ L inea r i z a t i on wi th 1 s t op po in t
∗∗∗∗∗∗∗∗∗∗∗∗∗ ’)
20 l i n=l i n e a r i z e ( ’Newmodel ’ , op po int (1 ) , i o ) ;
21 A=l i n . a ;
22 E( : , 1 )=eig (A) ;
23 A
24 % d i s p l a y ( ’ Abso lu te Eigen Values ’ )
25 % abs (E)
42
26 d i sp l ay ( ’ ∗∗∗∗∗∗∗∗∗∗ L i n e a r i z a t i o n with 2nd op point
∗∗∗∗∗∗∗∗∗∗∗∗∗ ’ )
27 l i n=l i n e a r i z e ( ’Newmodel ’ , op po int (2 ) , i o ) ;
28 A=l i n . a ;
29 E( : , 2 )=eig (A) ;
30 A
31 % d i s p l a y ( ’ Abso lu te Eigen Values ’ )
32 % abs (E)
33 d i sp l ay ( ’ ∗∗∗∗∗∗∗∗∗∗ L i n e a r i z a t i o n with 3 rd op po int
∗∗∗∗∗∗∗∗∗∗∗∗∗ ’ )
34 l i n=l i n e a r i z e ( ’Newmodel ’ , op po int (3 ) , i o ) ;
35 A=l i n . a ;
36 E( : , 3 )=eig (A) ;
37 A
38 % d i s p l a y ( ’ Abso lu te Eigen Values ’ )
39 % abs (E)
40 d i sp l ay ( ’ ∗∗∗∗∗∗∗∗∗∗ L i n e a r i z a t i o n with 4 th op po int
∗∗∗∗∗∗∗∗∗∗∗∗∗ ’ )
41 l i n=l i n e a r i z e ( ’Newmodel ’ , op po int (4 ) , i o ) ;
42 A=l i n . a ;
43 E( : , 4 )=eig (A) ;
44 A
45 % d i s p l a y ( ’ Abso lu te Eigen Values ’ )
46 % abs (E)
47 d i sp l ay ( ’ ∗∗∗∗∗∗∗∗∗∗ L i n e a r i z a t i o n with 5 th op po int
∗∗∗∗∗∗∗∗∗∗∗∗∗ ’ )
48 l i n=l i n e a r i z e ( ’Newmodel ’ , op po int (5 ) , i o ) ;
43
49 A=l i n . a ;
50 E( : , 5 )=eig (A) ;
51 A
52 % Set o f Eigen va l u e s
53 E
44
Code for control system characteristics
1 clc ;
2 %The cont inuous s t a t e−space model
3 A = [−0.4369 0 0 0 ;
4 0 .0159 −0.2390 0 0 ;
5 0 0 .0860 −0.1306 0 ;
6 0 .3880 0 .1200 0 .0860 −0.0330] ;
7 B= [−1 .554 ;
8 −0.032;
9 0 ;
10 1 . 5 8 6 ] ;
11 C= [ 1 0 0 0 ;
12 0 1 0 0 ;
13 0 0 1 0 ;
14 0 0 0 1 ] ;
15 Qc = ctrb (A,B) ; %c o n t r o l l a b i l i t y matrix
16 Qo = obsv (A,C) ; %o b s e r v a b i l i t y matrix
17 R1 = rank (Qc) ; %rank o f the c o n t r o l l a b i l i t y matrix
18 R2 = rank (Qo) ; %rank o f the o b s e r v a b i l i t y matrix
19 unob = length (A) − rank (Qo) ; %no . o f unobservab l e s t a t e s
20 uncon = length (A) − rank (Qc) ; %no . o f u n c on t r o l l a b l e s t a t e s
45
Code for linear quadratic regulation for first set of Q and R
1 x0 = [ 1 0 ; 8 ; 0 ; 0 ] ;
2 C=eye (4 ) ;
3 R = [ 3 0 8 0 0 ] ;
4 Q = [1000 50 50 50 ;
5 50 100 50 50 ;
6 50 50 50 50 ;
7 50 50 50 5 0 ] ;
8 F = 1∗eye (4 )
9 t0 =0;
10 t f =30;
11 tspan = [ t0 t f ] ;
12 A = [−0.4369 0 0 0 ;
13 0 .0159 −0.2390 0 0 ;
14 0 0 .0860 −0.1306 0 ;
15 0 .3880 0 .1200 0 .0860 −0.0330] ;
16 B = [−1 .554 ;
17 −0.032;
18 0 ;
19 1 . 5 8 6 ] ;
20 [ x , u , k]= l q r n s s (A,B, F ,Q,R, x0 , tspan ) ;
46
Code for linear quadratic regulation for second set of Q and R
1 x0 = [ 1 0 ; 8 ; 0 ; 0 ] ;
2 C=eye (4 ) ;
3 R = [ 5 0 0 0 ] ;
4 Q = [200 0 0 50 ;
5 0 50 50 50 ;
6 0 50 50 50 ;
7 50 50 50 5 0 ] ;
8 F = 1∗eye (4 )
9 t0 =0;
10 t f =30;
11 tspan = [ t0 t f ] ;
12 A = [−0.4369 0 0 0 ;
13 0 .0159 −0.2390 0 0 ;
14 0 0 .0860 −0.1306 0 ;
15 0 .3880 0 .1200 0 .0860 −0.0330] ;
16 B = [−1 .554 ;
17 −0.032;
18 0 ;
19 1 . 5 8 6 ] ;
20 [ x , u , k]= l q r n s s (A,B, F ,Q,R, x0 , tspan ) ;
47