# Optimal Tracking for a Divergent-Type Parabolic PDE System in Current Profile Control

Authors: S. Li; H. Jiang; Z. Ren; C. Xu
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

We consider the optimal tracking problem for a divergent-type parabolic PDE system, which can be used to model the spatial-temporal evolution of the magnetic diffusion process in a tokamak plasma. With the feedforward trajectories generated by numerical optimization, we derive a linear system to describe the error evolution. The weak variation approach is then extended to cylindrical coordinates to obtain a Riccati-type PDE for feedback kernel synthesis. Then, a finite difference method is used to give numerical solutions for the Riccati-type PDE. The obtained feedback kernel is finally used to simulate the closed-loop system. The simulation results demonstrate that the proposed tracking method can attenuate the tracking error effectively.


---

## Body

## 1. Introduction

Nuclear fusion, a process by which two nuclei fuse together to form a heavier nucleus, is accompanied by a release of energy, which is the result of the mass “lost” in the reaction. The amount of released energy is given by Einstein’s famous mass-energy formula. To make a fusion reaction possible, a certain amount of energy is required to bring two repellant nuclei carrying positive charges sufficiently close. To overcome the Coulomb barrier, the kinetic energy of the nuclei is increased by various heating methods (including ohmic heating, neutral beam injection (NBI), ion cyclotron resonance frequency (ICRF), electron cyclotron resonance frequency (ECRF), and lower hybrid heating (LHH)).The temperature required for a thermonuclear fusion reaction to take place is around 100 million degrees. At much lower temperatures (about 10 thousand degrees), the electrons and nuclei separate and create an ionized gas called plasma, also known as the fourth state of matter. The plasma conducts electricity and responds to magnetic fields, motivating a magnetic confinement approach to nuclear fusion. One type of magnetic confinement device is the tokamak (see Figure1), in which a torus-shaped intangible bottle is created by magnetic fields to confine the high-temperature plasma. It is well known that harvesting nuclear fusion is quite challenging. After more than fifty years of hard work, scientists now identify the challenges into three main aspects, that is, plasma physics, material science, and control (e.g., E. J. Synakowski’s presentation in the University Fusion Association (UFA) General Meeting, during the 51st American Physical Society Division of Plasma Physics (APS-DPP) Annual Meeting (November 2–6, 2009, Atlanta, Georgia, USA): Perspectives on Fusion Energy Research: On Our Science, Leverage and Credibility.) [1]. From the viewpoint of control systems, tokamak operation and control are very complex and strongly coupled with plasma dynamics and tokamak materials, such as magnetohydrodynamics stabilization, plasma confinement, steady-state operation, and purity control. A large amount of literature is now available at the interface of control systems and tokamak operation, such as the books [2, 3] and special issue papers (Refer to papers in the special issues titled “Control of Tokamak Plasmas: Part I” (October, 2005) and “Control of Tokamak Plasmas: Part II” (April, 2006) in IEEE Control Systems Magazine, organized by A. Pironti and M. Walker).Figure 1
Schematic of a tokamak chamber and magnetic profile.Feedback control of the main plasma macroscopic parameters (such as position, shape, total current, and line density) is reasonably well known [3]. However, control of internal profiles is still in its infancy [2]. Tokamaks are high order, distributed parameter, nonlinear systems with a large number of instabilities [4], so there are many extremely challenging mathematical modeling and control problems, which must be solved before a fusion power system becomes a viable entity [5–9]. The main focus of this paper is to study the current profile tracking of tokamak plasma which can be modeled by a parabolic PDE. In our previous work, we applied PDE-constrained optimization techniques to compute plasma discharge sequences minimizing a given cost function [7]. For example, the minimization problem can be a combination of the electricity consumption and practical-desired deviation at certain time instance.For the online implementation, we need to attenuate external perturbations and uncertainties. Therefore, feedback control is needed for the online implementation to guarantee that the system can track the preoptimized trajectories robustly as close as possible [10–12]. Thus, by defining a derivation dynamic system, we can formulate a feedback control problem of the derivation system which is governed by a linear parabolic PDE system. In this paper, we only use the boundary actuator for tracking due to the accuracy of magnetic actuation.Boundary control of parabolic PDEs has been investigated for a long time, mainly including the semigroup theory approach [13] and the recent infinite dimensional backstepping approach [14, 15]. The method used in this paper is motivated by [16] where a very neat and practical weak variation method is applied to derive the Riccati PDE for feedback kernel synthesis. Then, a routine spatial-temporal difference scheme is applied to solve the Riccati PDE. In this paper, we extend the method to solve a divergent-type parabolic PDE system model for the magnetic-flux profile which is a time-invariant system from rectangular coordinate system to polar coordinate system.We organize this paper as follows. In Section2, we state the optimal control problem for a divergent-type parabolic PDE model for the magnetic-flux profile with actuators at the boundary. In Section 3, we derive the optimal controller for the open-loop control PDE system using weak variation method. Further, we present the closed-loop optimal controller in Section 4. The numerical simulation results are presented in Section 5. Finally, Section 6 concludes the paper by stating the conclusions and further research topics.


## 2. Statement of the Optimal Control Problem

We consider a simplified divergent-type parabolic PDE model for the magnetic-flux profileψ:
(1)
∂
ψ
∂
t
(
r
,
t
)
=
1
r
∂
∂
r
(
r
D
(
r
,
t
)
∂
ψ
∂
r
(
r
,
t
)
)
+
η
(
r
,
t
)
j
(
r
,
t
)
,
with the Neumann boundary conditions
(2)
∂
ψ
∂
r
(
0
,
t
)
=
0
,
∂
ψ
∂
r
(
1
,
t
)
=
U
(
t
)
,
where ψ
(
r
,
t
) denotes the poloidal magnetic-flux around the tokamak; r denotes the normalized radius; t denotes time; U
(
t
) denotes the control variable acting at the right ending; D
(
r
,
t
) is divergent-type transport coefficient (which can be identified from the experimental data); η
(
r
,
t
) is the resistivity; and j
(
r
,
t
) is the noninductive current source. The initial condition for the magnetic flux is given by
(3)
ψ
(
r
,
0
)
=
ψ
0
(
r
)
.
Our goal is to design a state-feedback controller U
(
t
) that optimally regulates the system during the flat-top phase of the plasma discharge in tokamak plasmas. Thus, we wish to minimize the following cost functional over a finite time horizon:
(4)
min
⁡
U
⁡
J
(
ψ
,
U
)
≜
∫
0
T
‍
e
2
a
~
t
(
〈
ψ
(
r
,
t
)
-
ψ
des
⁡
(
r
,
t
)
,
Q
[
ψ
(
r
,
t
)
-
ψ
des
⁡
(
r
,
t
)
]
〉
≜
∫
0
T
‍
e
2
a
~
t
+
R
U
2
(
t
)
[
ψ
(
r
,
t
)
-
ψ
des
⁡
(
r
,
t
)
]
)
d
t
,
where a
~ is a positive constant, R is a positive weighting kernel, and
(5)
Q
ψ
(
r
,
t
)
≜
∫
0
1
‍
Q
(
r
,
y
)
ψ
(
y
,
t
)
y
d
y
,
(6)
〈
ψ
(
r
,
t
)
,
Q
ψ
(
r
,
t
)
〉
=
∫
0
1
‍
ψ
(
r
,
t
)
r
d
r
∫
0
1
‍
Q
(
r
,
y
)
ψ
(
y
,
t
)
y
d
y
.
Note that the desired profile ψ
des
⁡
(
r
,
t
) is generated by the following system:
(7)
∂
ψ
des
⁡
∂
t
(
r
,
t
)
=
1
r
∂
∂
r
(
r
D
(
r
,
t
)
∂
ψ
des
⁡
∂
r
(
r
,
t
)
)
+
η
(
r
,
t
)
j
¯
(
r
,
t
)
,
with the following Neumann boundary conditions:
(8)
∂
ψ
des
⁡
∂
r
(
0
,
t
)
=
0
,
∂
ψ
des
⁡
∂
r
(
1
,
t
)
=
U
¯
(
t
)
,
where j
¯
(
r
,
t
) and U
¯
(
t
) are the predefined internal and boundary inputs obtained by various numerical optimization approaches [7, 9, 12]. In addition, the open-loop optimal control input can be provided by combining the control parameterization method and gradient-based optimization (e.g., Loxton et al.’s paper [17]). The control problem is to synthesize a controller which can maintain the profile governed by (7)-(8). We assume that the internal input η
(
r
,
t
)
j
¯
(
r
,
t
) is fixed; that is, j
¯
(
r
,
t
)
=
j
(
r
,
t
). Then, we can define the new variables ψ
Δ
(
r
,
t
)
=
ψ
(
r
,
t
)
-
ψ
des
⁡
(
r
,
t
) and u
(
t
)
=
U
(
t
)
-
U
¯
(
t
), which are the state and control variables in the following optimal control problem:
(9)
min
⁡
u
⁡
J
(
ψ
Δ
(
r
,
t
)
,
u
(
t
)
)
=
∫
0
T
‍
e
2
a
~
t
(
〈
ψ
Δ
(
r
,
t
)
,
Q
ψ
Δ
(
r
,
t
)
〉
+
R
u
2
(
t
)
)
d
t
∂
ψ
Δ
∂
t
(
r
,
t
)
=
1
r
∂
∂
r
(
r
D
(
r
,
t
)
∂
ψ
Δ
∂
r
(
r
,
t
)
)
,
∂
ψ
Δ
∂
r
(
0
,
t
)
=
0
,
∂
ψ
Δ
∂
r
(
1
,
t
)
=
u
(
t
)
,
ψ
Δ
(
r
,
0
)
=
ψ
(
r
,
0
)
-
ψ
des
⁡
(
r
,
0
)
.
Without loss of generality, we can write ψ
Δ
(
r
,
t
) as ψ
(
r
,
t
). Then the optimal tracking problem becomes
(10)
min
⁡
u
⁡
J
(
ψ
(
r
,
t
)
,
u
(
t
)
)
=
∫
0
T
‍
e
2
a
~
t
(
〈
ψ
(
r
,
t
)
,
Q
ψ
(
r
,
t
)
〉
+
R
u
2
(
t
)
)
d
t
.
(11)
∂
ψ
∂
t
(
r
,
t
)
=
1
r
∂
∂
r
(
r
D
(
r
,
t
)
∂
ψ
∂
r
(
r
,
t
)
)
,
∂
ψ
∂
r
(
0
,
t
)
=
0
,
∂
ψ
∂
r
(
1
,
t
)
=
u
(
t
)
,
ψ
(
r
,
0
)
=
ψ
Δ
(
r
,
0
)
.
Now, we state our optimal control problem as follows.Problem
P
0. Given the PDE system (11), design an optimal state-feedback controller u
(
t
) such that the cost functional (10) is minimized.


## 3. Optimal Open-Loop Control

To solve ProblemP
0, we start by deriving the first order necessary conditions of the system for the open-loop finite-time problem. We have the following theorem.Theorem 1.
One considers the simplified poloidal magnetic-flux transport model (11). Let ψ
*
(
r
,
t
), u
*
(
t
), and λ
*
(
r
,
t
) denote the optimal state, control, and costate variables that minimize the cost functional (10). Then the optimality conditions are given as follows:
(12)
-
∂
λ
*
(
r
,
t
)
∂
t
=
D
(
r
,
t
)
∂
2
λ
*
(
r
,
t
)
∂
r
2
+
(
D
(
r
,
t
)
r
+
D
r
(
r
,
t
)
)
∂
λ
*
(
r
,
t
)
∂
r
+
2
e
2
a
~
t
Q
ψ
*
(
r
,
t
)
,
λ
r
*
(
0
,
t
)
=
0
,
λ
r
*
(
1
,
t
)
=
0
,
λ
*
(
r
,
T
)
=
0
,
u
*
(
t
)
=
-
1
2
R
e
-
2
a
~
t
D
(
1
,
t
)
λ
*
(
1
,
t
)
.Proof.
For simplicity, we will denote the optimal costate variablesλ
*
(
r
,
t
) by λ
(
r
,
t
) in the proof. The necessary optimality conditions are derived via weak variations motivated by Moura’s work [16]. Assume that the perturbed representations with respect to the optimal trajectories are given by
(13)
ψ
(
r
,
t
)
=
ψ
*
(
r
,
t
)
+
ϵ
δ
ψ
(
r
,
t
)
,
u
(
t
)
=
u
*
(
t
)
+
ϵ
δ
u
(
t
)
,
where δ
ψ
(
r
,
0
)
=
0, δ denotes the perturbation operator, and ϵ is an arbitrary constant. Then the perturbed cost function is
(14)
J
(
ψ
*
(
r
,
t
)
+
ϵ
δ
ψ
(
r
,
t
)
,
u
*
(
t
)
+
ϵ
δ
u
(
t
)
)
=
∫
0
T
‍
e
2
a
~
t
〈
ψ
*
(
r
,
t
)
+
ϵ
δ
ψ
(
r
,
t
)
,
h
h
h
h
h
h
h
h
h
h
h
Q
[
ψ
*
(
r
,
t
)
+
ϵ
δ
ψ
(
r
,
t
)
]
〉
d
t
+
∫
0
T
‍
e
2
a
~
t
R
(
u
*
(
t
)
+
ϵ
δ
u
(
t
)
)
2
d
t
.
We define a function G
(
ϵ
) of ϵ based on the perturbed cost functional (14) above and combine it with the dynamic system (11) by using the method of Lagrange multipliers as follows:
(15)
G
(
ϵ
)
=
J
(
ψ
*
(
r
,
t
)
+
ϵ
δ
ψ
(
r
,
t
)
,
u
*
(
t
)
+
ϵ
δ
u
(
t
)
)
+
∫
0
T
‍
〈
λ
(
r
,
t
)
,
D
(
r
,
t
)
∂
2
ψ
*
(
r
,
t
)
∂
r
2
h
h
h
h
h
h
+
ϵ
D
(
r
,
t
)
∂
2
δ
ψ
(
r
,
t
)
∂
r
2
〉
d
t
+
∫
0
T
〈
λ
(
r
,
t
)
,
(
D
r
(
r
,
t
)
+
D
(
r
,
t
)
r
)
∂
ψ
*
(
r
,
t
)
∂
r
h
h
h
h
h
h
+
ϵ
(
D
r
(
r
,
t
)
+
D
(
r
,
t
)
r
)
∂
δ
ψ
(
r
,
t
)
∂
r
〉
d
t
-
∫
0
T
‍
〈
λ
(
r
,
t
)
,
∂
∂
t
(
ψ
*
(
r
,
t
)
+
ϵ
δ
ψ
(
r
,
t
)
)
〉
d
t
,
where λ
(
r
,
t
) is the Lagrange multiplier. The necessary condition for optimality is
(16)
d
G
(
ϵ
)
d
ϵ
|
ϵ
=
0
=
0
.
A direct calculation gives
(17)
d
G
d
ϵ
(
ϵ
)
=
∫
0
T
‍
e
2
a
~
t
[
〈
δ
ψ
(
r
,
t
)
,
Q
[
ψ
*
(
r
,
t
)
+
ϵ
δ
ψ
(
r
,
t
)
]
〉
]
d
t
+
∫
0
T
‍
e
2
a
~
t
〈
ψ
*
(
r
,
t
)
+
ϵ
δ
ψ
(
r
,
t
)
,
Q
[
δ
ψ
(
r
,
t
)
]
〉
d
t
+
2
∫
0
T
‍
e
2
a
~
t
R
(
u
*
(
t
)
+
ϵ
δ
u
(
t
)
)
δ
u
(
t
)
d
t
+
∫
0
T
‍
〈
λ
(
r
,
t
)
,
D
(
r
,
t
)
∂
2
δ
ψ
(
r
,
t
)
∂
r
2
h
h
h
h
h
h
+
(
D
r
(
r
,
t
)
+
D
(
r
,
t
)
r
)
∂
δ
ψ
(
r
,
t
)
∂
r
∂
2
δ
ψ
(
r
,
t
)
∂
r
2
〉
d
t
-
∫
0
T
‍
〈
λ
(
r
,
t
)
,
∂
δ
ψ
(
r
,
t
)
∂
t
〉
d
t
.
By the definition of the inner product (6), we have
(18)
〈
λ
(
r
,
t
)
,
(
D
r
(
r
,
t
)
+
D
(
r
,
t
)
r
)
∂
δ
ψ
(
r
,
t
)
∂
r
〉
=
∫
0
1
‍
λ
(
r
,
t
)
(
D
r
(
r
,
t
)
r
+
D
(
r
,
t
)
)
∂
δ
ψ
(
r
,
t
)
∂
r
d
r
.
Using integration by parts, we can obtain
(19)
〈
λ
(
r
,
t
)
,
D
(
r
,
t
)
∂
2
δ
ψ
(
r
,
t
)
∂
r
2
〉
=
∫
0
1
‍
λ
(
r
,
t
)
D
(
r
,
t
)
∂
2
δ
ψ
(
r
,
t
)
∂
r
2
r
d
r
=
(
λ
(
r
,
t
)
D
(
r
,
t
)
r
∂
δ
ψ
(
r
,
t
)
∂
r
)
|
r
=
0
r
=
1
-
∫
0
1
‍
(
λ
(
r
,
t
)
D
(
r
,
t
)
+
λ
r
(
r
,
t
)
D
(
r
,
t
)
r
h
h
h
h
h
h
h
h
h
+
λ
(
r
,
t
)
D
r
(
r
,
t
)
r
)
×
∂
δ
ψ
(
r
,
t
)
∂
r
d
r
=
λ
(
1
,
t
)
D
(
1
,
t
)
∂
δ
ψ
(
r
,
t
)
∂
r
|
r
=
1
-
∫
0
1
‍
(
λ
(
r
,
t
)
D
(
r
,
t
)
+
λ
r
(
r
,
t
)
D
(
r
,
t
)
r
+
λ
D
r
(
r
,
t
)
r
)
h
h
h
h
h
h
×
∂
δ
ψ
(
r
,
t
)
∂
r
d
r
.
Note that ∂
δ
ψ
(
r
,
t
)
/
∂
r
|
r
=
1
=
δ
u
(
t
); hence by using (18) and (19), we obtain
(20)
〈
λ
(
r
,
t
)
,
D
(
r
,
t
)
∂
2
δ
ψ
(
r
,
t
)
∂
r
2
+
(
D
r
(
r
,
t
)
+
D
(
r
,
t
)
r
)
∂
δ
ψ
(
r
,
t
)
∂
r
∂
2
δ
ψ
(
r
,
t
)
∂
r
2
〉
=
λ
(
1
,
t
)
D
(
1
,
t
)
δ
u
(
t
)
-
∫
0
1
‍
λ
r
(
r
,
t
)
D
(
r
,
t
)
r
∂
δ
ψ
(
r
,
t
)
∂
r
d
r
=
λ
(
1
,
t
)
D
(
1
,
t
)
δ
u
(
t
)
-
λ
r
(
r
,
t
)
D
(
r
,
t
)
r
δ
ψ
(
r
,
t
)
|
r
=
0
r
=
1
+
∫
0
1
‍
(
λ
r
(
r
,
t
)
D
(
r
,
t
)
+
(
λ
r
r
(
r
,
t
)
D
(
r
,
t
)
+
λ
r
(
r
,
t
)
D
r
(
r
,
t
)
)
r
)
δ
ψ
(
r
,
t
)
d
r
=
λ
(
1
,
t
)
D
(
1
,
t
)
δ
u
(
t
)
-
λ
r
(
1
,
t
)
D
(
1
,
t
)
δ
ψ
(
1
,
t
)
+
〈
λ
r
(
r
,
t
)
D
(
r
,
t
)
r
+
λ
r
r
(
r
,
t
)
D
(
r
,
t
)
h
h
h
h
h
+
λ
r
(
r
,
t
)
D
r
(
r
,
t
)
,
δ
ψ
(
r
,
t
)
λ
r
(
r
,
t
)
D
(
r
,
t
)
r
〉
.
Moreover, we have
(21)
∫
0
T
‍
〈
λ
(
r
,
t
)
,
∂
δ
ψ
(
r
,
t
)
∂
t
〉
d
t
=
〈
λ
(
r
,
T
)
,
δ
ψ
(
r
,
T
)
〉
-
〈
λ
(
r
,
0
)
,
δ
ψ
(
r
,
0
)
〉
-
∫
0
T
‍
〈
∂
λ
(
r
,
t
)
∂
t
,
δ
ψ
(
r
,
t
)
〉
d
t
=
〈
λ
(
r
,
T
)
,
δ
ψ
(
r
,
T
)
〉
-
∫
0
T
‍
〈
∂
λ
(
r
,
t
)
∂
t
,
δ
ψ
(
r
,
t
)
〉
d
t
,
where δ
ψ
(
r
,
0
)
=
0 by definition. Substituting (19)–(21) into (17), we obtain
(22)
d
G
(
ϵ
)
d
ϵ
=
∫
0
T
‍
e
2
a
~
t
〈
δ
ψ
(
r
,
t
)
,
Q
[
(
ψ
*
(
r
,
t
)
+
ϵ
δ
ψ
(
r
,
t
)
)
]
〉
d
t
+
∫
0
T
‍
e
2
a
~
t
〈
ψ
*
(
r
,
t
)
+
ϵ
δ
ψ
(
r
,
t
)
,
Q
[
δ
ψ
(
r
,
t
)
]
〉
d
t
+
2
∫
0
T
‍
[
e
2
a
~
t
R
(
u
*
(
t
)
+
ϵ
δ
u
(
t
)
)
+
λ
(
1
,
t
)
D
(
1
,
t
)
]
×
δ
u
(
t
)
d
t
+
∫
0
T
‍
〈
λ
r
(
r
,
t
)
D
(
r
,
t
)
r
+
λ
r
r
(
r
,
t
)
D
(
r
,
t
)
h
h
h
h
h
h
+
λ
r
(
r
,
t
)
D
r
(
r
,
t
)
+
λ
t
,
δ
ψ
(
r
,
t
)
λ
r
(
r
,
t
)
D
(
r
,
t
)
r
〉
d
t
-
∫
0
T
‍
λ
r
(
1
,
t
)
D
(
1
,
t
)
δ
ψ
(
1
,
t
)
d
t
-
〈
λ
(
r
,
T
)
,
δ
ψ
(
r
,
T
)
〉
.
Thus,
(23)
d
G
(
ϵ
)
d
ϵ
|
ϵ
=
0
=
∫
0
T
‍
e
2
a
~
t
[
〈
δ
ψ
(
r
,
t
)
,
Q
[
ψ
*
(
r
,
t
)
]
〉
]
d
t
+
∫
0
T
‍
e
2
a
~
t
〈
ψ
*
(
r
,
t
)
,
Q
[
δ
ψ
(
r
,
t
)
]
〉
d
t
+
∫
0
T
‍
[
2
R
e
2
a
~
t
u
*
(
t
)
+
λ
(
1
,
t
)
D
(
1
,
t
)
]
δ
u
(
t
)
d
t
+
∫
0
T
‍
〈
λ
r
(
r
,
t
)
D
(
r
,
t
)
r
+
λ
r
r
(
r
,
t
)
D
(
r
,
t
)
h
h
h
h
h
h
+
λ
r
(
r
,
t
)
D
r
(
r
,
t
)
+
λ
t
,
δ
ψ
(
r
,
t
)
λ
r
(
r
,
t
)
D
(
r
,
t
)
r
〉
d
t
-
∫
0
T
‍
λ
r
(
1
,
t
)
D
(
1
,
t
)
δ
ψ
(
1
,
t
)
d
t
-
〈
λ
(
r
,
T
)
,
δ
ψ
(
r
,
T
)
〉
=
0
.
Note that we have used the following property:
(24)
〈
δ
ψ
(
r
,
t
)
,
Q
[
ψ
*
(
r
,
t
)
]
〉
=
〈
ψ
*
(
r
,
t
)
,
Q
[
δ
ψ
(
r
,
t
)
]
〉
.
Equality (23) holds for any δ
ψ
(
r
,
t
), δ
U
(
t
), and δ
ψ
(
r
,
T
) and therefore we get the following costate equation:
(25)
-
∂
λ
(
r
,
t
)
∂
t
=
D
(
r
,
t
)
∂
2
λ
(
r
,
t
)
∂
r
2
+
(
D
(
r
,
t
)
r
+
D
r
(
r
,
t
)
)
∂
λ
(
r
,
t
)
∂
r
+
2
e
2
a
~
t
Q
[
ψ
*
(
r
,
t
)
]
,
with the boundary conditions
(26)
λ
r
(
1
,
t
)
=
0
,
λ
r
(
0
,
t
)
=
0
,
and the final condition
(27)
λ
(
r
,
T
)
=
0
.
The optimal control input in open-loop form is given by
(28)
u
*
(
t
)
=
-
1
2
R
e
-
2
a
~
t
D
(
1
,
t
)
λ
(
1
,
t
)
.
This completes the proof.


## 4. Optimal State-Feedback Control

We now consider the closed-loop system with the control in state-feedback form. It is conventional in linear system theory to assume that the stateψ
(
r
,
t
) and the costate λ
(
r
,
t
) satisfy the following linear relationship [13]:
(29)
λ
(
r
,
t
)
=
P
ψ
(
r
,
t
)
=
∫
0
1
‍
P
(
r
,
y
,
t
)
ψ
(
y
,
t
)
y
d
y
,
where P
(
r
,
y
,
t
) is a function to be determined. We have the following theorem.Theorem 2.
The optimal control in state-feedback form is given by(30)
u
*
(
t
)
=
-
e
-
2
a
~
t
2
R
D
(
1
,
t
)
P
ψ
(
r
,
t
)
|
r
=
1
,
where P
(
r
,
y
,
t
) satisfies the following Riccati equation:
(31)
∂
P
(
r
,
y
,
t
)
∂
t
=
-
D
(
r
,
t
)
∂
2
P
(
r
,
y
,
t
)
∂
r
2
-
(
D
(
r
,
t
)
r
+
D
r
(
r
,
t
)
)
∂
P
(
r
,
y
,
t
)
∂
r
-
D
(
y
,
t
)
∂
2
P
(
r
,
y
,
t
)
∂
y
2
-
(
D
(
y
,
t
)
y
+
D
y
(
y
,
t
)
)
∂
P
(
r
,
y
,
t
)
∂
y
-
P
(
r
,
1
,
t
)
D
2
(
1
,
t
)
e
-
2
a
~
t
2
R
P
(
r
,
y
,
t
)
-
2
e
2
a
~
t
Q
(
r
,
y
,
t
)
,
with the boundary condition
(32)
P
r
(
0
,
y
,
t
)
=
P
r
(
1
,
y
,
t
)
=
P
y
(
r
,
0
,
t
)
=
P
y
(
r
,
1
,
t
)
=
0
and the final condition
(33)
P
(
r
,
y
,
T
)
=
0
.Proof.
By computing the time derivative ofλ
(
r
,
t
), we obtain
(34)
∂
λ
(
r
,
t
)
∂
t
=
∂
∂
t
∫
0
1
‍
P
(
r
,
y
,
t
)
ψ
(
y
,
t
)
y
d
y
=
∫
0
1
‍
[
P
t
(
r
,
y
,
t
)
ψ
(
y
,
t
)
+
P
(
r
,
y
,
t
)
ψ
t
(
y
,
t
)
]
y
d
y
.
Based on (11), we have
(35)
∫
0
1
‍
P
(
r
,
y
,
t
)
ψ
t
(
y
,
t
)
y
d
y
=
∫
0
1
‍
P
(
r
,
y
,
t
)
(
D
(
y
,
t
)
y
∂
2
ψ
(
y
,
t
)
∂
y
2
=
∫
0
1
‍
P
(
r
,
y
,
t
)
+
(
D
y
(
y
,
t
)
y
+
D
(
y
,
t
)
)
∂
ψ
(
y
,
t
)
∂
y
∂
2
ψ
(
y
,
t
)
∂
y
2
)
d
y
=
P
(
r
,
y
,
t
)
D
(
y
,
t
)
y
∂
ψ
(
y
,
t
)
∂
y
|
y
=
0
y
=
1
-
∫
0
1
‍
(
D
(
y
,
t
)
P
y
(
r
,
y
,
t
)
y
+
∂
D
(
y
,
t
)
∂
y
P
(
r
,
y
,
t
)
y
-
p
∫
0
1
‍
+
P
(
r
,
y
,
t
)
D
(
y
,
t
)
∂
D
(
y
,
t
)
∂
y
)
∂
ψ
(
y
,
t
)
∂
y
d
y
+
∫
0
1
‍
P
(
r
,
y
,
t
)
(
D
y
(
y
,
t
)
y
+
D
(
y
,
t
)
)
∂
ψ
(
y
,
t
)
∂
y
d
y
=
P
(
r
,
1
,
t
)
D
(
1
,
t
)
∂
ψ
(
1
,
t
)
∂
y
-
∫
0
1
‍
D
(
y
,
t
)
y
P
y
(
r
,
y
,
t
)
∂
ψ
(
y
,
t
)
∂
y
d
y
=
P
(
r
,
1
,
t
)
D
(
1
,
t
)
∂
ψ
(
1
,
t
)
∂
y
-
P
y
(
r
,
y
,
t
)
D
(
y
,
t
)
y
ψ
(
y
,
t
)
|
y
=
0
y
=
1
+
∫
0
1
‍
(
P
y
y
(
r
,
y
,
t
)
D
(
y
,
t
)
y
+
P
y
(
r
,
y
,
t
)
D
y
(
y
,
t
)
y
+
P
y
(
r
,
y
,
t
)
D
(
y
,
t
)
)
ψ
(
y
,
t
)
d
y
=
P
(
r
,
1
,
t
)
D
(
1
,
t
)
u
(
t
)
-
P
y
(
r
,
1
,
t
)
D
(
1
,
t
)
ψ
(
1
,
t
)
+
∫
0
1
‍
(
P
y
y
(
r
,
y
,
t
)
D
(
y
,
t
)
y
+
P
y
(
r
,
y
,
t
)
D
y
(
y
,
t
)
y
+
P
y
(
r
,
y
,
t
)
D
(
y
,
t
)
)
ψ
(
y
,
t
)
d
y
.
We choose the following boundary condition for P
(
r
,
y
,
t
):
(36)
P
(
r
,
0
,
t
)
=
0
,
P
y
(
r
,
1
,
t
)
=
0
,
P
y
(
r
,
0
,
t
)
=
0
.
We then immediately obtain
(37)
∫
0
1
‍
P
(
r
,
y
,
t
)
ψ
t
(
y
,
t
)
y
d
y
=
P
(
r
,
1
,
t
)
D
(
1
,
t
)
u
(
t
)
+
∫
0
1
‍
(
P
y
y
(
r
,
y
,
t
)
D
(
y
,
t
)
y
+
P
y
(
r
,
y
,
t
)
D
y
(
y
,
t
)
y
h
h
h
h
h
h
h
h
h
+
P
y
(
r
,
y
,
t
)
D
(
y
,
t
)
)
ψ
(
y
,
t
)
d
y
=
P
(
r
,
1
,
t
)
D
2
(
1
,
t
)
e
-
2
a
~
t
2
R
λ
(
1
,
t
)
+
∫
0
1
‍
(
P
y
y
(
r
,
y
,
t
)
D
(
y
,
t
)
y
+
P
y
(
r
,
y
,
t
)
D
y
(
y
,
t
)
y
h
h
h
h
h
h
h
h
h
+
P
y
(
r
,
y
,
t
)
D
(
y
,
t
)
)
ψ
(
y
,
t
)
d
y
=
P
(
r
,
1
,
t
)
D
2
(
1
,
t
)
e
-
2
a
~
t
2
R
(
∫
0
1
‍
P
(
r
,
y
,
t
)
ψ
(
y
,
t
)
y
d
y
)
+
∫
0
1
‍
(
P
y
y
(
r
,
y
,
t
)
D
(
y
,
t
)
y
+
P
y
(
r
,
y
,
t
)
D
y
(
y
,
t
)
y
h
h
h
h
h
h
h
h
h
+
P
y
(
r
,
y
,
t
)
D
(
y
,
t
)
)
ψ
(
y
,
t
)
d
y
.
Therefore, we have
(38)
∂
λ
(
r
,
t
)
∂
t
=
∫
0
1
‍
(
P
t
(
r
,
y
,
t
)
ψ
(
y
,
t
)
+
P
(
r
,
y
,
t
)
ψ
t
(
y
,
t
)
)
y
d
y
=
∫
0
1
‍
P
t
(
r
,
y
,
t
)
ψ
(
y
,
t
)
y
d
y
+
P
(
r
,
1
,
t
)
D
2
(
1
,
t
)
e
-
2
a
~
t
2
R
(
∫
0
1
‍
P
(
r
,
y
,
t
)
ψ
(
y
,
t
)
y
d
y
)
+
∫
0
1
‍
(
P
y
y
(
r
,
y
,
t
)
D
(
y
,
t
)
y
+
P
y
(
r
,
y
,
t
)
D
y
(
y
,
t
)
y
h
h
h
h
h
h
h
h
h
+
P
y
(
r
,
y
,
t
)
D
(
y
,
t
)
)
ψ
(
y
,
t
)
d
y
.
Recalling optimality conditions (12) in Theorem 1, we obtain
(39)
∂
λ
(
r
,
t
)
∂
t
=
-
D
(
r
,
t
)
∫
0
1
‍
∂
2
P
(
r
,
y
,
t
)
∂
r
2
ψ
(
y
,
t
)
y
d
y
-
(
D
(
r
,
t
)
r
+
D
r
(
r
,
t
)
)
∫
0
1
‍
∂
P
(
r
,
y
,
t
)
∂
r
ψ
(
y
,
t
)
y
d
y
-
2
e
2
a
~
t
∫
0
1
‍
Q
(
r
,
y
,
t
)
ψ
(
y
,
t
)
y
d
y
.
Combining (38) and (39), then we obtain
(40)
∂
P
(
r
,
y
,
t
)
∂
t
=
-
D
(
r
,
t
)
∂
2
P
(
r
,
y
,
t
)
∂
r
2
-
(
D
(
r
,
t
)
r
+
D
r
(
r
,
t
)
)
∂
P
(
r
,
y
,
t
)
∂
r
-
D
(
y
,
t
)
∂
2
P
(
r
,
y
,
t
)
∂
y
2
-
(
D
(
y
,
t
)
y
+
D
y
(
y
,
t
)
)
∂
P
(
r
,
y
,
t
)
∂
y
-
P
(
r
,
1
,
t
)
D
2
(
1
,
t
)
e
-
2
a
~
t
2
R
P
(
r
,
y
,
t
)
-
2
e
2
a
~
t
Q
(
r
,
y
)
.
By virtue of the optimality conditions given in (12), λ
(
r
,
t
) must satisfy λ
r
(
0
,
t
)
=
0, λ
r
(
1
,
t
)
=
0, and λ
(
r
,
T
)
=
0 and hence
(41)
λ
r
(
0
,
t
)
=
∫
0
1
‍
P
r
(
0
,
y
,
t
)
ψ
(
y
,
t
)
y
d
y
=
0
,
λ
r
(
1
,
t
)
=
∫
0
1
‍
P
r
(
1
,
y
,
t
)
ψ
(
y
,
t
)
y
d
y
=
0
,
λ
(
r
,
T
)
=
∫
0
1
‍
P
(
r
,
y
,
T
)
ψ
(
y
,
T
)
y
d
y
=
0
,
which hold for any ψ
(
y
,
t
) and thus we obtain the following boundary conditions:
(42)
P
r
(
0
,
y
,
t
)
=
0
,
P
r
(
1
,
y
,
t
)
=
0
,
ϕ
r
(
0
,
t
)
=
0
,
ϕ
r
(
1
,
t
)
=
0
,
and the final condition
(43)
P
(
r
,
y
,
T
)
=
0
.
This completes the proof.


## 5. Numerical Simulations

In this section, we simulate the results to verify the effectiveness of the proposed method. The coefficient functions and parameters for the system are given in Table1. Our numerical simulation based on the finite difference method is carried out within the MATLAB programming environment.Table 1
Coefficient functions and parameters for numerical simulation example.
Transport coefficient
D
(
r
,
t
)
=
e
r
t
Initial condition
ψ
(
r
,
0
)
=
sin
⁡
⁡
π
r
Terminal time
T
=
2First, we examine the open-loop system described by (11); that is, u
(
t
)
=
0. The evolution of the open-loop system without feedback control is shown in Figure 2. Second, the designed optimal feedback controller (30) is applied in order to check the closed-loop system (11); see Figure 3. The simulation results show that, with the designed optimal feedback controller, the original state is suppressed and the state converges to zero which verifies the effectiveness of our proposed method. The boundary control input is shown in Figure 4.Figure 2
Open-loop responseψ
(
r
,
t
).Figure 3
Closed-loop responseψ
(
x
,
t
) corresponding to the optimal controller u
*
(
t
) (30).Figure 4
Boundary control input forψ
(
1
,
t
).We only present simulation results in terms of the parameters and functions in Table1 but this method will work for other cases with parameters and functions satisfying conditions in Theorem 2.


## 6. Conclusion and Future Work

In this paper, an infinite dimensional feedback controller is obtained for the optimal tracking problem arising in the flat-top phase of plasma discharges in a tokamak. Our method avoids discretizing the process model into lumped parameter systems (i.e., ODEs) but directly uses infinite dimensional LQ control method to synthesize a feedback controller. To provide the feedback kernel numerically, the finite difference method is used to solve the Riccati-type PDE. Although only simulation results have been presented in the current work, our future work will apply this approach to trajectory tracking in the EAST (the Experimental Advanced Superconducting Tokamak, Institute of Plasma Physics, Chinese Academy of Sciences, Hefei, Anhui, China) device.



---

*Source: 940965-2014-06-11.xml*