# Stability Analysis of Pneumatic Cabin Pressure Regulating System with Complex Nonlinear Characteristics

Authors: Xinhua Zheng; Lili Xie; Lizhuo Liu
Journal: Journal of Control Science and Engineering
 (2015)

---

## Abstract

Stability of pneumatic cabin pressure regulating system with complex nonlinear characteristics is considered. The mathematical model of each component is obtained and given in detail. The governing equations of the considered system consist of 8 differential equations. In the circumstance, commonly used methods of nonlinear system analysis are not applicable. Therefore a new method is proposed to construct phase plane trajectories numerically. The calculation steps are given in detail. And convergence region of numerical calculation and limits on step size is defined. The method is applied constructing phase plane trajectories for considered cabin pressure regulating system. Phase plane analysis shows that there exists a limit cycle, which is responsible for pressure pulsating in aircraft cabin. After parameters adjustment, excellent stability characteristics are acquired. And the validity of this method is confirmed by the simulation.


---

## Body

## 1. Introduction

It is known that cabin pressure regulating system is used to provide a safe and comfortable surrounding for crew and passengers at high altitude [1]. Now cabin pressure regulating system is necessary for commercial airplanes and military aircrafts cruising at high altitude. For commercial airplanes, digital pressurization systems had become a new developing trend. And in some newly designed airplanes, it is integrated with other systems [2]. The implement of this new technology advantages improvement in intellectualization and automation of control. At the same time, it improves comfort of the passenger and reduces labor cost in maintenance work. But these numerical systems are electrical energy dependent and sensitive to electromagnetic interference. Therefore in military aircrafts, pneumatic cabin pressure regulating systems are still preferable.The cabin pressure regulating system considered in present paper is an all pneumatic system. It regulates cabin pressure automatically and no human interference is needed. It works pneumatically and needs no electrical energy. This makes it a system of high reliability, even if electromagnetic interference of high intensive cannot have any influence on its performance. And it adapts automatically for use on high altitude airports. It also has additional desirable characteristics of light weight, small volume, and low maintenance. These characteristics satisfy exactly the needs of advanced military airplanes. Now it is used widely in newly developed military airplanes in China.But at certain flying conditions, pulsating of cabin pressure arises. The pulsating of cabin pressure is undesirable for it may lead to common symptoms of airplane ear including hearing loss, discomfort in the ear, or mild to moderate pain. Continuous existence of such symptoms distracts attention of pilots and reduces their work efficiency or even leads to permanent hearing impairment [3]. Therefore the stability of cabin pressure must be ensured. The main object of the present paper is to define the cause of the cabin pressure pulsating and eliminate it.


## 2. Cabin Pressure Regulating System and Its Dynamic Equation

### 2.1. Diagram of Cabin Pressure Regulating System

As Figure1 displays, the considered cabin pressure regulating system comprises five main parts, schedule pressure generator, rate limiter, pneumatic amplifier, outflow valve, and cabin. Comprised of two independent components, schedule pressure generator is the core of the system. Its main function is to produce cabin pressure schedule signal to control outflow valve. The rate limiter limits maximum rising rate and maximum dropping rate of scheduled pressure. The output signal of rate limiter is too weak to drive pneumatic outflow valve; therefore pneumatic amplifier is used to amplify its driven capability. And as a pneumatically sourced element, it is capable of driving outflow valve directly. Outflow valve is the actuator of the system. Its function is to open valve in accordance with pressure in its control chamber, which is a pneumatic load of amplifier. At certain altitude, the flow rate of cabin air to ambient circumstance is directly proportional to the opening of outflow valve. And by this means cabin pressure is regulated.Figure 1
The diagram of cabin pressure regulating system.


### 2.2. Mathematical Model of Components

The differential equations of the cabin pressure regulating system are built on basis of the following assumptions [4].(1)
The cabin temperatureT
p keeps invariant, and T
p = 20°C.
(2)
The cabin volume keeps constant.
(3)
The air filling the cabin is ideal gas.
(4)
The leakage area of the cabin is constant.Because the cabin temperature keeps constant, the flow quantity can be measured in volume in the mathematical model built for each component. This greatly simplifies computation.In the following mathematical models,p
1, p
2, p
3, and p
4 are the output pressures of schedule pressure generator, rate limiter, pneumatic amplifier, and cabin, l
1 is the opening of valve 1 in the active component of schedule pressure generator, l
31, l
32 are the opening of valve 31 and valve 32, which connect pneumatic amplifier to higher and lower pressure source, respectively, l
4 is the opening of outflow valve, s
1, s
31, s
32, and s
4 are the corresponding effective flow area of valve 1, valve 31, valve 32, and outflow valve, and P
a is the air pressure around the aircraft.


#### 2.2.1. Mathematical Model for Schedule Pressure Generator

Schedule pressure generator is composed of two components: proportional pressure generator and differential pressure generator. The two components take part in schedule pressure modulating at different altitude of flight. Proportional pressure generator works at lower altitude of flight, and differential pressure generator takes over its job a height about 5100 meter. The two components have similar structure and similar governing equations. Therefore in the present paper, only differential pressure generator is analyzed. Correspondingly only stability at altitude above 5100 meter is considered. The stability at lower altitude can be analyzed similarly.The functional structure of differential pressure generator is shown in Figure2. Having a double nozzle-flapper structure disposed in modulating chamber, differential pressure generator has a similar structure as proportional pressure generator [5]. The differential equation of its valve’s opening l
1 is
(1)
l
˙
1
=
(
p
a
-
p
1
c
)
A
1
+
F
10
+
k
1
l
1
f
1
,
0
<
l
1
<
L
M
1
,
0
,
l
1
≤
0
or
l
1
≥
L
M
1
,
where A
1 is the effective area of pressure applied on diaphragms, F
10 is the preload spring force on diaphragms, k
1 is the equivalent elasticity coefficient of the springs, and f
1 is the kinetic friction coefficient of valve 1 movement.Figure 2
Functional structure of differential pressure generator.In the situation of dynamic balance, volumes flowing in and out of feedback chamber are equal. According to flow rate calculation formula equation (3) in [6] and equation (1) in [7], p
1
c is dynamic solution of the following equation:
(2)
q
c
1
p
1
c
1
-
p
3
/
p
1
c
-
b
1
1
-
b
1
2
=
q
c
2
p
4
1
-
p
1
c
/
p
4
-
b
2
1
-
b
2
2
,
where q
c
1 is the saturation volume flow of capillary pipe mc1, q
c
2 is the saturation volume flow of capillary pipe mc2, b
1 is the critical pressure ratio of capillary mc1, and b
2 is the critical pressure ratio of capillary mc2. Therefore p
1
c is a variable which lies between p
3 and p
4.In accordance with (15) in [8], the change rate of p
1 is
(3)
p
˙
1
=
240
s
11
p
4
p
1
p
4
-
p
1
V
1
-
240
S
10
-
s
11
p
1
p
3
-
C
1
p
1
-
p
3
+
C
1
V
1
,
where s
11
=
c
·
l
1 is the effective flow area of valve 1 and S
10 is the maximum effective flow area of valve 1. C
1 is a constant, which stands for the pressure drop of P
3. V
1 is the volume of modulating chamber of differential pressure generator.


#### 2.2.2. Mathematical Model for Rate Limiter

The rate limiter is composed of a capillary pipe mc3 connecting to a gas container. By limits volume flowing through capillary pipe, the change rate of output pressure in gas container is limited. Dynamic equation ofp
2 is as follows:
(4)
p
˙
2
=
q
c
3
p
1
V
2
1
-
p
2
/
p
1
-
b
3
1
-
b
3
2
,
p
2
<
p
1
,
-
q
c
3
p
2
V
2
1
-
p
2
/
p
1
-
b
3
1
-
b
3
2
,
p
2
>
p
1
,
where q
c
3 is the saturation volume flow of capillary pipe mc3, b
3 is the critical pressure ratio of capillary pipe mc3, and V
2 is the volume of gas container.


#### 2.2.3. Mathematical Model for Pneumatic Amplifier

The functional structure of pneumatic amplifier is shown in Figure3. It also has a similar structure as proportional pressure generator. In pneumatic amplifier, there are two independent valves: valve 31 and valve 32. Suppose the opening of valve 31 and valve 32 are l
31 and l
32, respectively, the relevant differential equations are as follows:
(5)
l
˙
31
=
(
p
2
-
p
3
)
A
3
-
F
310
-
k
31
l
31
f
31
,
0
<
l
31
<
L
M
31
,
0
,
l
31
≤
0
or
l
31
≥
L
M
31
,
l
˙
32
=
(
p
3
-
p
2
)
A
3
-
F
320
-
k
32
l
32
f
32
,
0
<
l
32
<
L
M
32
,
0
,
l
32
≤
0
or
l
32
≥
L
M
32
,
p
˙
3
=
240
s
31
p
4
p
3
p
4
-
p
3
V
3
-
120
s
32
p
3
V
3
,
p
a
p
3
≤
0.5
,
240
s
31
p
4
p
3
p
4
-
p
3
V
3
-
240
s
32
p
3
p
a
p
3
-
p
a
V
3
,
p
a
p
3
>
0.5
,
where s
31
=
c
·
l
31 is the dynamic charging area of valve 31, s
32
=
c
·
l
32 is the dynamic discharging area of valve 32, and V
3 is the total volume of modulating chamber and pneumatic load of pneumatic amplifier.Figure 3
Functional structure of pneumatic amplifier.


#### 2.2.4. Mathematical Model for Outflow Valve and Cabin Pressure

Cabin pressure and its change rate is the control object of cabin pressure regulating system. The velocity of the outflow valve opening is(6)
l
˙
4
=
(
p
4
-
p
3
)
A
4
-
F
40
-
k
4
l
4
f
4
,
0
<
l
4
<
L
M
4
,
0
,
l
4
≤
0
or
l
4
≥
L
M
4
,
where A
4 is the effective area of air pressure applied on diaphragm of outflow valve, F
40 is the force preload on outflow valve by reset spring in control chamber of outflow valve, k
4 is the elasticity coefficient of the reset spring, and f
4 is the kinetic friction coefficient of outflow valve movement.The volume flow rate of outflow valveG
Out is
(7)
G
Out
=
120
s
4
p
4
,
p
a
p
4
≤
0.5
240
s
4
p
4
p
a
p
4
1
-
p
a
p
4
,
p
a
p
4
>
0.5
.According to (13) in [9] and (2) in [10], the change rate of cabin pressure is
(8)
p
˙
4
=
G
S
(
h
,
h
˙
)
-
G
L
-
G
Out
V
C
·
p
4
,
where V
C is the volume of aircraft cabin, G
L is the volume flow rate caused by cabin leakage, and G
S
(
h
,
h
˙
) is the volume flow into the cabin caused by air supply system.


## 2.1. Diagram of Cabin Pressure Regulating System

As Figure1 displays, the considered cabin pressure regulating system comprises five main parts, schedule pressure generator, rate limiter, pneumatic amplifier, outflow valve, and cabin. Comprised of two independent components, schedule pressure generator is the core of the system. Its main function is to produce cabin pressure schedule signal to control outflow valve. The rate limiter limits maximum rising rate and maximum dropping rate of scheduled pressure. The output signal of rate limiter is too weak to drive pneumatic outflow valve; therefore pneumatic amplifier is used to amplify its driven capability. And as a pneumatically sourced element, it is capable of driving outflow valve directly. Outflow valve is the actuator of the system. Its function is to open valve in accordance with pressure in its control chamber, which is a pneumatic load of amplifier. At certain altitude, the flow rate of cabin air to ambient circumstance is directly proportional to the opening of outflow valve. And by this means cabin pressure is regulated.Figure 1
The diagram of cabin pressure regulating system.


## 2.2. Mathematical Model of Components

The differential equations of the cabin pressure regulating system are built on basis of the following assumptions [4].(1)
The cabin temperatureT
p keeps invariant, and T
p = 20°C.
(2)
The cabin volume keeps constant.
(3)
The air filling the cabin is ideal gas.
(4)
The leakage area of the cabin is constant.Because the cabin temperature keeps constant, the flow quantity can be measured in volume in the mathematical model built for each component. This greatly simplifies computation.In the following mathematical models,p
1, p
2, p
3, and p
4 are the output pressures of schedule pressure generator, rate limiter, pneumatic amplifier, and cabin, l
1 is the opening of valve 1 in the active component of schedule pressure generator, l
31, l
32 are the opening of valve 31 and valve 32, which connect pneumatic amplifier to higher and lower pressure source, respectively, l
4 is the opening of outflow valve, s
1, s
31, s
32, and s
4 are the corresponding effective flow area of valve 1, valve 31, valve 32, and outflow valve, and P
a is the air pressure around the aircraft.


### 2.2.1. Mathematical Model for Schedule Pressure Generator

Schedule pressure generator is composed of two components: proportional pressure generator and differential pressure generator. The two components take part in schedule pressure modulating at different altitude of flight. Proportional pressure generator works at lower altitude of flight, and differential pressure generator takes over its job a height about 5100 meter. The two components have similar structure and similar governing equations. Therefore in the present paper, only differential pressure generator is analyzed. Correspondingly only stability at altitude above 5100 meter is considered. The stability at lower altitude can be analyzed similarly.The functional structure of differential pressure generator is shown in Figure2. Having a double nozzle-flapper structure disposed in modulating chamber, differential pressure generator has a similar structure as proportional pressure generator [5]. The differential equation of its valve’s opening l
1 is
(1)
l
˙
1
=
(
p
a
-
p
1
c
)
A
1
+
F
10
+
k
1
l
1
f
1
,
0
<
l
1
<
L
M
1
,
0
,
l
1
≤
0
or
l
1
≥
L
M
1
,
where A
1 is the effective area of pressure applied on diaphragms, F
10 is the preload spring force on diaphragms, k
1 is the equivalent elasticity coefficient of the springs, and f
1 is the kinetic friction coefficient of valve 1 movement.Figure 2
Functional structure of differential pressure generator.In the situation of dynamic balance, volumes flowing in and out of feedback chamber are equal. According to flow rate calculation formula equation (3) in [6] and equation (1) in [7], p
1
c is dynamic solution of the following equation:
(2)
q
c
1
p
1
c
1
-
p
3
/
p
1
c
-
b
1
1
-
b
1
2
=
q
c
2
p
4
1
-
p
1
c
/
p
4
-
b
2
1
-
b
2
2
,
where q
c
1 is the saturation volume flow of capillary pipe mc1, q
c
2 is the saturation volume flow of capillary pipe mc2, b
1 is the critical pressure ratio of capillary mc1, and b
2 is the critical pressure ratio of capillary mc2. Therefore p
1
c is a variable which lies between p
3 and p
4.In accordance with (15) in [8], the change rate of p
1 is
(3)
p
˙
1
=
240
s
11
p
4
p
1
p
4
-
p
1
V
1
-
240
S
10
-
s
11
p
1
p
3
-
C
1
p
1
-
p
3
+
C
1
V
1
,
where s
11
=
c
·
l
1 is the effective flow area of valve 1 and S
10 is the maximum effective flow area of valve 1. C
1 is a constant, which stands for the pressure drop of P
3. V
1 is the volume of modulating chamber of differential pressure generator.


### 2.2.2. Mathematical Model for Rate Limiter

The rate limiter is composed of a capillary pipe mc3 connecting to a gas container. By limits volume flowing through capillary pipe, the change rate of output pressure in gas container is limited. Dynamic equation ofp
2 is as follows:
(4)
p
˙
2
=
q
c
3
p
1
V
2
1
-
p
2
/
p
1
-
b
3
1
-
b
3
2
,
p
2
<
p
1
,
-
q
c
3
p
2
V
2
1
-
p
2
/
p
1
-
b
3
1
-
b
3
2
,
p
2
>
p
1
,
where q
c
3 is the saturation volume flow of capillary pipe mc3, b
3 is the critical pressure ratio of capillary pipe mc3, and V
2 is the volume of gas container.


### 2.2.3. Mathematical Model for Pneumatic Amplifier

The functional structure of pneumatic amplifier is shown in Figure3. It also has a similar structure as proportional pressure generator. In pneumatic amplifier, there are two independent valves: valve 31 and valve 32. Suppose the opening of valve 31 and valve 32 are l
31 and l
32, respectively, the relevant differential equations are as follows:
(5)
l
˙
31
=
(
p
2
-
p
3
)
A
3
-
F
310
-
k
31
l
31
f
31
,
0
<
l
31
<
L
M
31
,
0
,
l
31
≤
0
or
l
31
≥
L
M
31
,
l
˙
32
=
(
p
3
-
p
2
)
A
3
-
F
320
-
k
32
l
32
f
32
,
0
<
l
32
<
L
M
32
,
0
,
l
32
≤
0
or
l
32
≥
L
M
32
,
p
˙
3
=
240
s
31
p
4
p
3
p
4
-
p
3
V
3
-
120
s
32
p
3
V
3
,
p
a
p
3
≤
0.5
,
240
s
31
p
4
p
3
p
4
-
p
3
V
3
-
240
s
32
p
3
p
a
p
3
-
p
a
V
3
,
p
a
p
3
>
0.5
,
where s
31
=
c
·
l
31 is the dynamic charging area of valve 31, s
32
=
c
·
l
32 is the dynamic discharging area of valve 32, and V
3 is the total volume of modulating chamber and pneumatic load of pneumatic amplifier.Figure 3
Functional structure of pneumatic amplifier.


### 2.2.4. Mathematical Model for Outflow Valve and Cabin Pressure

Cabin pressure and its change rate is the control object of cabin pressure regulating system. The velocity of the outflow valve opening is(6)
l
˙
4
=
(
p
4
-
p
3
)
A
4
-
F
40
-
k
4
l
4
f
4
,
0
<
l
4
<
L
M
4
,
0
,
l
4
≤
0
or
l
4
≥
L
M
4
,
where A
4 is the effective area of air pressure applied on diaphragm of outflow valve, F
40 is the force preload on outflow valve by reset spring in control chamber of outflow valve, k
4 is the elasticity coefficient of the reset spring, and f
4 is the kinetic friction coefficient of outflow valve movement.The volume flow rate of outflow valveG
Out is
(7)
G
Out
=
120
s
4
p
4
,
p
a
p
4
≤
0.5
240
s
4
p
4
p
a
p
4
1
-
p
a
p
4
,
p
a
p
4
>
0.5
.According to (13) in [9] and (2) in [10], the change rate of cabin pressure is
(8)
p
˙
4
=
G
S
(
h
,
h
˙
)
-
G
L
-
G
Out
V
C
·
p
4
,
where V
C is the volume of aircraft cabin, G
L is the volume flow rate caused by cabin leakage, and G
S
(
h
,
h
˙
) is the volume flow into the cabin caused by air supply system.


## 2.2.1. Mathematical Model for Schedule Pressure Generator

Schedule pressure generator is composed of two components: proportional pressure generator and differential pressure generator. The two components take part in schedule pressure modulating at different altitude of flight. Proportional pressure generator works at lower altitude of flight, and differential pressure generator takes over its job a height about 5100 meter. The two components have similar structure and similar governing equations. Therefore in the present paper, only differential pressure generator is analyzed. Correspondingly only stability at altitude above 5100 meter is considered. The stability at lower altitude can be analyzed similarly.The functional structure of differential pressure generator is shown in Figure2. Having a double nozzle-flapper structure disposed in modulating chamber, differential pressure generator has a similar structure as proportional pressure generator [5]. The differential equation of its valve’s opening l
1 is
(1)
l
˙
1
=
(
p
a
-
p
1
c
)
A
1
+
F
10
+
k
1
l
1
f
1
,
0
<
l
1
<
L
M
1
,
0
,
l
1
≤
0
or
l
1
≥
L
M
1
,
where A
1 is the effective area of pressure applied on diaphragms, F
10 is the preload spring force on diaphragms, k
1 is the equivalent elasticity coefficient of the springs, and f
1 is the kinetic friction coefficient of valve 1 movement.Figure 2
Functional structure of differential pressure generator.In the situation of dynamic balance, volumes flowing in and out of feedback chamber are equal. According to flow rate calculation formula equation (3) in [6] and equation (1) in [7], p
1
c is dynamic solution of the following equation:
(2)
q
c
1
p
1
c
1
-
p
3
/
p
1
c
-
b
1
1
-
b
1
2
=
q
c
2
p
4
1
-
p
1
c
/
p
4
-
b
2
1
-
b
2
2
,
where q
c
1 is the saturation volume flow of capillary pipe mc1, q
c
2 is the saturation volume flow of capillary pipe mc2, b
1 is the critical pressure ratio of capillary mc1, and b
2 is the critical pressure ratio of capillary mc2. Therefore p
1
c is a variable which lies between p
3 and p
4.In accordance with (15) in [8], the change rate of p
1 is
(3)
p
˙
1
=
240
s
11
p
4
p
1
p
4
-
p
1
V
1
-
240
S
10
-
s
11
p
1
p
3
-
C
1
p
1
-
p
3
+
C
1
V
1
,
where s
11
=
c
·
l
1 is the effective flow area of valve 1 and S
10 is the maximum effective flow area of valve 1. C
1 is a constant, which stands for the pressure drop of P
3. V
1 is the volume of modulating chamber of differential pressure generator.


## 2.2.2. Mathematical Model for Rate Limiter

The rate limiter is composed of a capillary pipe mc3 connecting to a gas container. By limits volume flowing through capillary pipe, the change rate of output pressure in gas container is limited. Dynamic equation ofp
2 is as follows:
(4)
p
˙
2
=
q
c
3
p
1
V
2
1
-
p
2
/
p
1
-
b
3
1
-
b
3
2
,
p
2
<
p
1
,
-
q
c
3
p
2
V
2
1
-
p
2
/
p
1
-
b
3
1
-
b
3
2
,
p
2
>
p
1
,
where q
c
3 is the saturation volume flow of capillary pipe mc3, b
3 is the critical pressure ratio of capillary pipe mc3, and V
2 is the volume of gas container.


## 2.2.3. Mathematical Model for Pneumatic Amplifier

The functional structure of pneumatic amplifier is shown in Figure3. It also has a similar structure as proportional pressure generator. In pneumatic amplifier, there are two independent valves: valve 31 and valve 32. Suppose the opening of valve 31 and valve 32 are l
31 and l
32, respectively, the relevant differential equations are as follows:
(5)
l
˙
31
=
(
p
2
-
p
3
)
A
3
-
F
310
-
k
31
l
31
f
31
,
0
<
l
31
<
L
M
31
,
0
,
l
31
≤
0
or
l
31
≥
L
M
31
,
l
˙
32
=
(
p
3
-
p
2
)
A
3
-
F
320
-
k
32
l
32
f
32
,
0
<
l
32
<
L
M
32
,
0
,
l
32
≤
0
or
l
32
≥
L
M
32
,
p
˙
3
=
240
s
31
p
4
p
3
p
4
-
p
3
V
3
-
120
s
32
p
3
V
3
,
p
a
p
3
≤
0.5
,
240
s
31
p
4
p
3
p
4
-
p
3
V
3
-
240
s
32
p
3
p
a
p
3
-
p
a
V
3
,
p
a
p
3
>
0.5
,
where s
31
=
c
·
l
31 is the dynamic charging area of valve 31, s
32
=
c
·
l
32 is the dynamic discharging area of valve 32, and V
3 is the total volume of modulating chamber and pneumatic load of pneumatic amplifier.Figure 3
Functional structure of pneumatic amplifier.


## 2.2.4. Mathematical Model for Outflow Valve and Cabin Pressure

Cabin pressure and its change rate is the control object of cabin pressure regulating system. The velocity of the outflow valve opening is(6)
l
˙
4
=
(
p
4
-
p
3
)
A
4
-
F
40
-
k
4
l
4
f
4
,
0
<
l
4
<
L
M
4
,
0
,
l
4
≤
0
or
l
4
≥
L
M
4
,
where A
4 is the effective area of air pressure applied on diaphragm of outflow valve, F
40 is the force preload on outflow valve by reset spring in control chamber of outflow valve, k
4 is the elasticity coefficient of the reset spring, and f
4 is the kinetic friction coefficient of outflow valve movement.The volume flow rate of outflow valveG
Out is
(7)
G
Out
=
120
s
4
p
4
,
p
a
p
4
≤
0.5
240
s
4
p
4
p
a
p
4
1
-
p
a
p
4
,
p
a
p
4
>
0.5
.According to (13) in [9] and (2) in [10], the change rate of cabin pressure is
(8)
p
˙
4
=
G
S
(
h
,
h
˙
)
-
G
L
-
G
Out
V
C
·
p
4
,
where V
C is the volume of aircraft cabin, G
L is the volume flow rate caused by cabin leakage, and G
S
(
h
,
h
˙
) is the volume flow into the cabin caused by air supply system.


## 3. Scheme of Stability Analysis and Its Realization

### 3.1. Expression of the System in General Form

Sety
1
=
p
1, y
2
=
p
2, y
3
=
p
3, y
4
=
p
4, y
5
=
l
1, y
6
=
l
31, y
7
=
l
32, and y
8
=
l
4, and according to (1)~(8), write their differential expressions in a general form [11], as shown in
(9)
d
y
1
d
t
=
f
1
(
t
,
y
1
,
y
2
,
…
,
y
8
)
d
y
2
d
t
=
f
2
(
t
,
y
1
,
y
2
,
…
,
y
8
)
⋮
d
y
8
d
t
=
f
8
t
,
y
1
,
y
2
,
…
,
y
8
,
and if the initial condition is given in
(10)
y
1
0
=
α
1
,
y
2
0
=
α
2
,
…
,
y
8
0
=
α
8
,
then stability analysis of the system described by differential equation set (9) with initial condition (10) is the object of the present paper.


### 3.2. Equilibrium Point of the System

An equilibrium point of system (9) with initial condition (10) is a solution, y
*
=
(
y
1
*
,
y
2
*
,
…
,
y
8
*
), of the following algebraic equation:
(11)
f
1
(
t
,
y
1
*
,
y
2
*
,
…
,
y
8
*
)
=
0
f
2
(
t
,
y
1
*
,
y
2
*
,
…
,
y
8
*
)
=
0
⋮
f
8
(
t
,
y
1
*
,
y
2
*
,
…
,
y
8
*
)
=
0
.At certain height (ht), the considered cabin pressure regulating system has only one equilibrium point. Suppose ht is higher than switching height, at which differential pressure component takes over the modulating of schedule pressure. The equilibrium point of the system (9) at height ht is
(12)
l
1
*
=
1
2
L
Max
⁡
l
2
*
=
l
3
*
=
0
,
p
1
*
=
p
2
*
=
p
3
*
=
p
ht
+
Δ
s
d
=
p
a
+
Δ
s
d
,
p
4
*
=
20
F
10
+
p
ht
+
0.5
k
1
L
M
1
s
-
19
p
3
*
=
20
F
10
+
p
a
+
0.5
k
1
L
M
1
s
-
19
p
3
*
,
l
4
*
=
(
p
4
*
-
p
3
*
)
A
4
-
F
40
k
4
,
where p
ht is pressure of atmosphere at height ht, p
ht
=
p
a.


### 3.3. Definition of Stability

The stability in present paper is a kind of global uniform asymptotic stability, and its definition [12] is as follows.System is said to beasymptotically stable about its equilibrium if it is stable in the sense of Lyapunov, and furthermore, there exists a constant, δ
=
δ
(
t
0
)
>
0, such that
(13)
x
(
t
0
)
<
δ
⟹
x
t
⟶
0
,
as
t
⟶
∞
.
Theasymptotic stability is said to beuniform if the existing constant is independent of t
0 over [
0
,
∞
] and is said to beglobal if the convergence x
(
t
)
→
0 is independent of the initial state x over the entire domain on which the system is defined.On the basis of the above definition, the cabin pressure regulating system has its special requirement on stability; that is, the transient must complete within 25 seconds.


### 3.4. Scheme of Stability Analysis and Its Realization

As stated above, the considered cabin pressure control system is a nonlinear system consisting of 8 differential equations and form of each equation is very complicated. Therefore stability analysis of the system is difficult.Describing function method, phase plane method, and Lyapunov method are the three methods commonly used to analyze stability of systems. But to a complex system containing multiple nonlinearities, describing function method is very complicated and unique for special problem [13]. Therefore it is not a universal method. Lyapunov method is an important technique that is used for not only analyzing but also designing stability of the system. And it can be applied conveniently to systems containing new control techniques and systems with special characteristics [14, 15]. But the difficulty of using Lyapunov function lies in founding and screening Lyapunov function. In this circumstance, Lyapunov function methods cannot be used as it is impossible to find a V function. And classic meaning of phase plane method, designed for second order systems, cannot solve the problem too. Therefore in the present paper, a numerical method is proposed and applied constructing phase plane trajectories by approximating time-domain solution of the system. Then on this basis, system stability can be analyzed in accordance with phase plane method.For System (9) with initial condition (10), within certain time interval 0
<
t
<
T, there are 8 functions y
1
t
,
y
2
t
,
…
,
y
8
(
t
) that satisfy each of the differential equations together with all the initial conditions. But it is difficult to find such functions analytically. Therefore we should first approximate y
1
t
,
y
2
t
,
…
,
y
8
(
t
) numerically at specified time intervals. The existence and uniqueness of the time domain solution are assured by Definition 5.16 and Theorem 5.17 in [11].To approximate time-domain solution of a complex nonlinear system subject to initial condition, the uniqueness of time-domain solution needs to be ensured firstly. Then a suitable numerical algorithm should be selected and adapted for special needs of considered problem. And selected step size should fulfill requirements of numerical stability of the calculation.Therefore the secondary object of present paper is to propose a new method to construct phase plane numerically, which extends the scope of phase plane analysis from two-dimension to multidimension.If the system of first-order differential equations (9), subject to the initial conditions (10), has a unique solution y
1
t
,
y
2
t
,
…
,
y
m
(
t
), then numerical methods could be used approximating its time-domain solution within certain period of time.Taking into consideration computational complexity and accuracy, the improved Euler method is selected, and its computational formulas are as follows:(14)
y
i
(
n
+
1
)
=
y
i
n
+
h
2
k
i
1
+
k
i
2
,
k
i
1
=
f
t
n
,
y
i
n
,
k
i
2
=
f
t
n
+
h
,
y
i
n
+
h
k
i
1
,
for
each
n
=
0
,
1
,
…
,
N
,
i
=
1
,
2
,
…
,
8
,
where N is number of steps needed to completed computation at time interval [
0
,
T
].In order to construct phase plane trajectories for considered cabin pressure regulating system, approximated time domain solution as well as its derivative needs to be acquired. Therefore the computation formula is adapted, as shown in(15)
y
i
(
n
+
1
)
=
y
i
n
+
h
2
k
i
1
+
k
i
2
,
k
i
1
=
f
t
n
,
y
i
n
,
k
i
2
=
f
t
n
+
h
,
y
i
n
+
h
k
i
1
,
y
˙
i
(
n
+
1
)
=
f
t
n
+
1
,
y
i
n
+
1
,
for
each
n
=
0
,
1
,
…
,
N
,
i
=
1
,
2
,
…
,
8
.In present paper an m-file is programmed on MATLAB platform to acquire data needed, according to equation set (15). Use the notation M
y
i
j and D
y
i
j for each j
=
0,1
,
…
,
N and i
=
1,2
,
…
,
8, to denote an approximation to y
i
(
t
j
) and its derivative y
i
′
(
t
j
), respectively. Set
(16)
M
y
11
=
α
1
,
M
y
21
=
α
2
,
…
,
M
y
81
=
α
8
.
The process of computation is described as follows. (1)
Set number of cycleN in accordance with step size h and time limit T, N
=
T
/
h.
(2)
Set initial value as(17)
M
y
1,1
=
α
1
,
M
y
2,1
=
α
2
,
…
,
M
y
8,1
=
α
8
,
D
y
1,1
=
f
0
,
M
y
1,1
,
D
y
2,1
=
f
0
,
M
y
2,1
,
…
,
D
y
8,1
=
f
0
,
M
y
8,1
.
(3)
Whiles increase from 1 to N, carry out computation from Step 4 to Step 7:
(4)
(18)
k
11
k
21
⋮
k
81
=
D
y
1
s
D
y
2
s
⋮
D
y
8
s
,
(5)
(19)
k
12
k
22
⋮
k
82
=
f
1
t
s
+
h
,
M
y
s
+
h
k
11
f
2
t
s
+
h
,
M
y
s
+
h
k
21
⋮
f
8
t
s
+
h
,
M
y
s
+
h
k
81
,
(6)
(20)
M
y
1
(
s
+
1
)
M
y
2
(
s
+
1
)
⋮
M
y
8
(
s
+
1
)
=
M
y
1
s
M
y
2
s
⋮
M
y
8
s
+
h
2
k
11
+
k
12
k
21
+
k
22
⋮
k
81
+
k
82
,
(7)
(21)
D
y
1
(
s
+
1
)
D
y
2
(
s
+
1
)
⋮
D
y
8
(
s
+
1
)
=
f
1
t
s
+
1
,
M
y
1
(
s
+
1
)
f
2
t
s
+
1
,
M
y
2
(
s
+
1
)
⋮
f
8
t
s
+
1
,
M
y
8
(
s
+
1
)
.
(8)
Plot interestedith pair of variable (M
y
i
,
D
y
i) of acquired data (here i may be any number between 1 and 8). Phase plane trajectories are now constructed, and stability of system with initial-value may be analyzed in accordance with phase plane methods.


### 3.5. Convergence of Calculation

Define model equation of absolute stability as linear differential equation:(22)
y
˙
1
y
˙
2
⋮
y
˙
8
=
∂
f
1
∂
y
1
0
⋯
0
0
∂
f
2
∂
y
2
⋯
0
⋮
⋮
⋱
⋮
0
0
…
∂
f
8
∂
y
8
y
1
y
2
⋮
y
8
=
λ
1
0
⋯
0
0
λ
2
⋯
0
⋮
⋮
⋱
⋮
0
0
⋯
λ
8
y
1
y
2
⋮
y
8
.The region of absolute stability for the improved Euler method is1
+
h
λ
+
(
h
λ
)
2
/
2
<
1, and the interval of absolute stability is -
2
<
h
λ
<
0, for each of λ
1
,
λ
2
,
…
,
λ
8.Setλ
i
X
=
max
⁡
0
≤
t
≤
T
∂
f
i
/
∂
y
i, λ
X
=
max
⁡
(
λ
i
X
), (
i
=
1,2
,
…
,
8
); then to ensure the results of calculation approximate to the time-domain solution y
(
t
), the interval of  h  value is limited by 0
<
h
<
2
/
λ
X.


## 3.1. Expression of the System in General Form

Sety
1
=
p
1, y
2
=
p
2, y
3
=
p
3, y
4
=
p
4, y
5
=
l
1, y
6
=
l
31, y
7
=
l
32, and y
8
=
l
4, and according to (1)~(8), write their differential expressions in a general form [11], as shown in
(9)
d
y
1
d
t
=
f
1
(
t
,
y
1
,
y
2
,
…
,
y
8
)
d
y
2
d
t
=
f
2
(
t
,
y
1
,
y
2
,
…
,
y
8
)
⋮
d
y
8
d
t
=
f
8
t
,
y
1
,
y
2
,
…
,
y
8
,
and if the initial condition is given in
(10)
y
1
0
=
α
1
,
y
2
0
=
α
2
,
…
,
y
8
0
=
α
8
,
then stability analysis of the system described by differential equation set (9) with initial condition (10) is the object of the present paper.


## 3.2. Equilibrium Point of the System

An equilibrium point of system (9) with initial condition (10) is a solution, y
*
=
(
y
1
*
,
y
2
*
,
…
,
y
8
*
), of the following algebraic equation:
(11)
f
1
(
t
,
y
1
*
,
y
2
*
,
…
,
y
8
*
)
=
0
f
2
(
t
,
y
1
*
,
y
2
*
,
…
,
y
8
*
)
=
0
⋮
f
8
(
t
,
y
1
*
,
y
2
*
,
…
,
y
8
*
)
=
0
.At certain height (ht), the considered cabin pressure regulating system has only one equilibrium point. Suppose ht is higher than switching height, at which differential pressure component takes over the modulating of schedule pressure. The equilibrium point of the system (9) at height ht is
(12)
l
1
*
=
1
2
L
Max
⁡
l
2
*
=
l
3
*
=
0
,
p
1
*
=
p
2
*
=
p
3
*
=
p
ht
+
Δ
s
d
=
p
a
+
Δ
s
d
,
p
4
*
=
20
F
10
+
p
ht
+
0.5
k
1
L
M
1
s
-
19
p
3
*
=
20
F
10
+
p
a
+
0.5
k
1
L
M
1
s
-
19
p
3
*
,
l
4
*
=
(
p
4
*
-
p
3
*
)
A
4
-
F
40
k
4
,
where p
ht is pressure of atmosphere at height ht, p
ht
=
p
a.


## 3.3. Definition of Stability

The stability in present paper is a kind of global uniform asymptotic stability, and its definition [12] is as follows.System is said to beasymptotically stable about its equilibrium if it is stable in the sense of Lyapunov, and furthermore, there exists a constant, δ
=
δ
(
t
0
)
>
0, such that
(13)
x
(
t
0
)
<
δ
⟹
x
t
⟶
0
,
as
t
⟶
∞
.
Theasymptotic stability is said to beuniform if the existing constant is independent of t
0 over [
0
,
∞
] and is said to beglobal if the convergence x
(
t
)
→
0 is independent of the initial state x over the entire domain on which the system is defined.On the basis of the above definition, the cabin pressure regulating system has its special requirement on stability; that is, the transient must complete within 25 seconds.


## 3.4. Scheme of Stability Analysis and Its Realization

As stated above, the considered cabin pressure control system is a nonlinear system consisting of 8 differential equations and form of each equation is very complicated. Therefore stability analysis of the system is difficult.Describing function method, phase plane method, and Lyapunov method are the three methods commonly used to analyze stability of systems. But to a complex system containing multiple nonlinearities, describing function method is very complicated and unique for special problem [13]. Therefore it is not a universal method. Lyapunov method is an important technique that is used for not only analyzing but also designing stability of the system. And it can be applied conveniently to systems containing new control techniques and systems with special characteristics [14, 15]. But the difficulty of using Lyapunov function lies in founding and screening Lyapunov function. In this circumstance, Lyapunov function methods cannot be used as it is impossible to find a V function. And classic meaning of phase plane method, designed for second order systems, cannot solve the problem too. Therefore in the present paper, a numerical method is proposed and applied constructing phase plane trajectories by approximating time-domain solution of the system. Then on this basis, system stability can be analyzed in accordance with phase plane method.For System (9) with initial condition (10), within certain time interval 0
<
t
<
T, there are 8 functions y
1
t
,
y
2
t
,
…
,
y
8
(
t
) that satisfy each of the differential equations together with all the initial conditions. But it is difficult to find such functions analytically. Therefore we should first approximate y
1
t
,
y
2
t
,
…
,
y
8
(
t
) numerically at specified time intervals. The existence and uniqueness of the time domain solution are assured by Definition 5.16 and Theorem 5.17 in [11].To approximate time-domain solution of a complex nonlinear system subject to initial condition, the uniqueness of time-domain solution needs to be ensured firstly. Then a suitable numerical algorithm should be selected and adapted for special needs of considered problem. And selected step size should fulfill requirements of numerical stability of the calculation.Therefore the secondary object of present paper is to propose a new method to construct phase plane numerically, which extends the scope of phase plane analysis from two-dimension to multidimension.If the system of first-order differential equations (9), subject to the initial conditions (10), has a unique solution y
1
t
,
y
2
t
,
…
,
y
m
(
t
), then numerical methods could be used approximating its time-domain solution within certain period of time.Taking into consideration computational complexity and accuracy, the improved Euler method is selected, and its computational formulas are as follows:(14)
y
i
(
n
+
1
)
=
y
i
n
+
h
2
k
i
1
+
k
i
2
,
k
i
1
=
f
t
n
,
y
i
n
,
k
i
2
=
f
t
n
+
h
,
y
i
n
+
h
k
i
1
,
for
each
n
=
0
,
1
,
…
,
N
,
i
=
1
,
2
,
…
,
8
,
where N is number of steps needed to completed computation at time interval [
0
,
T
].In order to construct phase plane trajectories for considered cabin pressure regulating system, approximated time domain solution as well as its derivative needs to be acquired. Therefore the computation formula is adapted, as shown in(15)
y
i
(
n
+
1
)
=
y
i
n
+
h
2
k
i
1
+
k
i
2
,
k
i
1
=
f
t
n
,
y
i
n
,
k
i
2
=
f
t
n
+
h
,
y
i
n
+
h
k
i
1
,
y
˙
i
(
n
+
1
)
=
f
t
n
+
1
,
y
i
n
+
1
,
for
each
n
=
0
,
1
,
…
,
N
,
i
=
1
,
2
,
…
,
8
.In present paper an m-file is programmed on MATLAB platform to acquire data needed, according to equation set (15). Use the notation M
y
i
j and D
y
i
j for each j
=
0,1
,
…
,
N and i
=
1,2
,
…
,
8, to denote an approximation to y
i
(
t
j
) and its derivative y
i
′
(
t
j
), respectively. Set
(16)
M
y
11
=
α
1
,
M
y
21
=
α
2
,
…
,
M
y
81
=
α
8
.
The process of computation is described as follows. (1)
Set number of cycleN in accordance with step size h and time limit T, N
=
T
/
h.
(2)
Set initial value as(17)
M
y
1,1
=
α
1
,
M
y
2,1
=
α
2
,
…
,
M
y
8,1
=
α
8
,
D
y
1,1
=
f
0
,
M
y
1,1
,
D
y
2,1
=
f
0
,
M
y
2,1
,
…
,
D
y
8,1
=
f
0
,
M
y
8,1
.
(3)
Whiles increase from 1 to N, carry out computation from Step 4 to Step 7:
(4)
(18)
k
11
k
21
⋮
k
81
=
D
y
1
s
D
y
2
s
⋮
D
y
8
s
,
(5)
(19)
k
12
k
22
⋮
k
82
=
f
1
t
s
+
h
,
M
y
s
+
h
k
11
f
2
t
s
+
h
,
M
y
s
+
h
k
21
⋮
f
8
t
s
+
h
,
M
y
s
+
h
k
81
,
(6)
(20)
M
y
1
(
s
+
1
)
M
y
2
(
s
+
1
)
⋮
M
y
8
(
s
+
1
)
=
M
y
1
s
M
y
2
s
⋮
M
y
8
s
+
h
2
k
11
+
k
12
k
21
+
k
22
⋮
k
81
+
k
82
,
(7)
(21)
D
y
1
(
s
+
1
)
D
y
2
(
s
+
1
)
⋮
D
y
8
(
s
+
1
)
=
f
1
t
s
+
1
,
M
y
1
(
s
+
1
)
f
2
t
s
+
1
,
M
y
2
(
s
+
1
)
⋮
f
8
t
s
+
1
,
M
y
8
(
s
+
1
)
.
(8)
Plot interestedith pair of variable (M
y
i
,
D
y
i) of acquired data (here i may be any number between 1 and 8). Phase plane trajectories are now constructed, and stability of system with initial-value may be analyzed in accordance with phase plane methods.


## 3.5. Convergence of Calculation

Define model equation of absolute stability as linear differential equation:(22)
y
˙
1
y
˙
2
⋮
y
˙
8
=
∂
f
1
∂
y
1
0
⋯
0
0
∂
f
2
∂
y
2
⋯
0
⋮
⋮
⋱
⋮
0
0
…
∂
f
8
∂
y
8
y
1
y
2
⋮
y
8
=
λ
1
0
⋯
0
0
λ
2
⋯
0
⋮
⋮
⋱
⋮
0
0
⋯
λ
8
y
1
y
2
⋮
y
8
.The region of absolute stability for the improved Euler method is1
+
h
λ
+
(
h
λ
)
2
/
2
<
1, and the interval of absolute stability is -
2
<
h
λ
<
0, for each of λ
1
,
λ
2
,
…
,
λ
8.Setλ
i
X
=
max
⁡
0
≤
t
≤
T
∂
f
i
/
∂
y
i, λ
X
=
max
⁡
(
λ
i
X
), (
i
=
1,2
,
…
,
8
); then to ensure the results of calculation approximate to the time-domain solution y
(
t
), the interval of  h  value is limited by 0
<
h
<
2
/
λ
X.


## 4. Results

The method stated above is applied to construct phase plane trajectories of cabin pressure regulating system with 2 initial conditions.At heightht
=
7000 m, the first initial condition is given as
(23)
p
10
=
80
,
p
20
=
80
,
p
30
=
80
,
p
40
=
80.5
,
l
10
=
0
,
l
310
=
0
,
l
320
=
0
,
l
40
=
0
.And at heightht
=
7000 m, the second initial condition is
(24)
p
10
=
78
,
p
20
=
78
,
p
30
=
78
,
p
40
=
79.5
,
l
10
=
0
,
l
310
=
0
,
l
320
=
0
,
l
40
=
0
.According to Section3.3, the calculated λ
X
<
308.4; then the interval of h value is 0
<
h
<
0.0064. Then step size h is chosen as h
=
0.001s. Take into consideration the actual requirement of the system, set T
=
30s, and then N
=
T
/
h
=
30000. And data is acquired by an m-file programmed in accordance with computation process listed in Section 3.2.In the cabin pressure regulating system, pressure stability ofp
2 and p
4 is variables interested. The phase plane trajectories of p
2 and p
4 are shown in Figures 4 and 5, respectively.Figure 4
Phase plane trajectories ofp
2.Figure 5
Phase plane trajectories ofp
4.Figure4 shows that, under the 2 given initial conditions, a limit cycle exists in pressure p
2. That means pressure p
2 will not be stable. The fluctuation in p
2 will pass on and be reproduced in cabin pressure, may be in reduced or comparable amplitude. Figure 5 shows an amplitude reduced limit cycle exists in p
4 phase plane trajectories. This is not allowed for the fluctuation may be aggravated in the complex environment of the flying aircrafts.In order to ensure stability of cabin pressurep
4, the parameters of the system should be adjusted. According to design experience, enhancing energy storage capability in the pressure regulating process helps to realize stability of system. And reducing dead zone threshold of pneumatic amplifier helps in eliminating limit cycle. Therefore, by substituting new type of polymer material for old one, and using new type of spring with bigger elastic coefficient in differential pressure generator, some system parameters are changed. And preload forces on amplifier’s valve are reduced slightly. New set of parameters is as follows:
(25)
k
1
=
4
N
/
cm
,
k
31
=
0.2
N
/
cm
,
k
32
=
0.2
N
/
cm
,
f
1
=
0.1
,
f
31
=
f
32
=
0.1
,
f
4
=
0.5
,
F
1
=
9.4
N
,
F
310
=
0.15
N
,
F
320
=
0.15
N
,
F
4
=
1.5
N
.And the corresponding replaced old set of parameters is as follows:(26)
k
1
=
3.5
N
/
cm
,
k
31
=
0.17
N
/
cm
,
k
32
=
0.17
N
/
cm
,
f
1
=
0.06
,
f
31
=
f
32
=
0.05
,
f
4
=
0.5
,
F
1
=
9.4
N
,
F
310
=
0.2
N
,
F
320
=
0.2
N
,
F
4
=
1.5
N
.After parameters adjustment, the new phase plane trajectories ofp
4 are shown in Figure 6. Figure 6 shows that no limit cycle exists in the system and the phase plane trajectories goes to the same equilibrium point even though the initial conditions are different. The simulating result shown in Figure 7 confirms this conclusion. Figure 7 shows that the stability of cabin pressure is achieved quickly and little fluctuation will happen. On this base, we come to the conclusion that the system is stable after parameter modification.Figure 6
Phase plane trajectories of cabin pressurep
4 after parameters modification.Figure 7
Simulation result of cabin pressurep
4.It needs to point out that, comparing with two-dimensional system, in a multidimensional system some difference exists in their phase plane trajectories; that is, the trajectories may come across limit cycle before it finally enters into it, as shown by Figure4.


## 5. Conclusion

The considered cabin pressure regulating system in present paper is a multidimensional complex nonlinear system. To analyze stability of cabin pressure regulating system which is used in advanced military airplane, a numerical method of construct phase plane trajectories is proposed in the present paper. This method extends the scope of phase plane analysis from two-dimension to multidimension. The method is applied to stability analysis of cabin pressure regulating system, and the validity of analysis is verified by simulation result. The method is applicable to complex nonlinear system of any dimension and may provide an effective way for stability analysis of complex nonlinear systems. It is instructive for researchers facing similar problems.



---

*Source: 871621-2015-02-24.xml*