Finite rank distributed control for the resistive diffusion equation using damping assignment

A first extension of the IDA-PBC control synthesis to infinite dimensional port Hamiltonian systems is investigated, using the same idea as for the finite dimensional case, that is transform the original model into a closed loop target Hamiltonian model using feedback control. To achieve this goal both finite rank distributed control and boundary control are used. The proposed class of considered port Hamiltonian distributed parameters systems is first defined. Then the matching equation is derived for this class before considering the particular case of damping assignment on the resistive diffusion example, for the radial diffusion of the poloidal magnetic flux in tokamak reactors. 1. A tribute to Professor Abdelhaq El Jai. The merging of Mediterranean cultures gave to the humanity thoroughly generous scholars and great humanist thinkers. Professor Abdelhaq El jai is certainly a brilliant representative of this long lineage. We would like here to pay a tribute to his outstanding contributions as a scientific, a humanist leader, a community builder and a cheerful friend. Certainly his original ideas on distributed parameters systems and control opened the way for numerous research and participated greatly to create a living scientific community, noticeably fed with his deep insights into the spatial (and regional) properties of these systems. Nevertheless, these ideas would never have created such an enthusiastic and friendly community without his exceptional human qualities. Personally, among all these qualities, we were greatly impressed by Professor El Jai kindness, benevolence, rigor and self-requirement, tolerance, enthusiasm and confidence to others’s capacities. He is a living example for all of us about how to lead our career and live our life as honorable women and men. 2. Introduction. Infinite dimensional port-Hamiltonian systems has recently become more and more popular either in the systems theory community as a class of naturally well-posed (linear) systems [3] or suitable for control designs based on Casimir functionals [11], Control by Interconnection (CbI) [8] or energy shaping [6]. On the other hand, Interconnection and Damping Assignment Passivity Based Control (IDA-PBC) methods have been successful in the control of nonlinear (and linear) finite dimensional Port-Controlled Hamiltonian (PCH) systems [9, 10]. It has also been used extensively in the control of finite dimensional PCH systems obtained via geometric spatial reduction of port-Hamiltonian systems, such as in the works on plasma current control for tokamaks developed previously by the authors [19, 20]. Roughly speaking, it makes use of the feedback control to match the

original system with a desired system written in the form of a closed loop asymptotically stable Hamiltonian systems. In this paper we propose a first extension of this idea to infinite dimensional open port-Hamiltonian systems. More precisely, while keeping the geometrical interconnection structure (namely the Dirac structure) unchanged, we propose energy shaping and damping assignment to match a restricted class of closed loop port-Hamiltonian systems. In order to achieve this result, we need to use the finite rank distributed control in the state equation and the boundary control simultaneously.
Indeed our application case concerns the control of the radial profile of the magnetic flux in tokamak toroidal reactors. In this example both boundary (coils) and distributed (antennas) actuators exist. The non-inductive current injection plays the role of distributed control while the loop voltage creates a boundary action. The distributed control is finite rank since only the total incoming power and the angle for the injected waves are controlled while the radial distribution (shape) of the control action is fixed for a given actuator. Unlike in the "traditional" boundary control methods where the boundary action is homogenized in system state evolution equation, the feedback distributed control is used to match a system with a homogeneous state equation and a new boundary control action which includes both the original control and the propagation of the distributed control.
The paper is organized as follows. In section 2 we define the class of controlled port-Hamiltonian systems, the class of target systems and the resulting matching equation for the control design. In section 3 we apply this idea to the example of the resistive diffusion equation for the poloidal flux in tokamaks. Both distributed non inductive current and boundary loop voltage are used as control variables. Damping assignment is performed in order to achieve some prescribed dissipation and the resulting asymptotic stability. Finally the feasibility of this control is discussed through some simulation examples and a robustness analysis.
3.1. The class of considered original and target port-Hamiltonian systems. This paper investigates the control problem for a class of distributed parameters port Hamiltonian systems with both boundary and distributed finite rank controls defined as: where the distributed control u 1 is just a time dependent scalar signal u 1 (t). The control constraint g (x) represents the spatial distribution of this control action. J = −J * is a formally skew-adjoint differential operator (cf. [7] Corollary 3.3). For the sake of simplicity we will consider spatial operators of the form J = P 1 ∂ z + P 0 , where P 1 ∈ M n (R) is non-singular symmetric matrix and where P 0 = −P T 0 ∈ M n (R) is a skew-symmetric one, although this class may be generalized to higher order spatial derivatives [5,14]. Denoting Z = [a, b], the spatial domain, the state space or space of configurations is chosen as x ∈ X = L 2 ([0, +∞) × Z). The dissipation is defined using the non negative self-adjoint operator R 0, and the total energy stored in the system using the hamiltonian smooth function H : where H is the energy density. B is a differential operator induced on the boundary ∂Z = {a, b} by the differential operator J in the sense of the following lemma 3.1 in [7], which is revisited as follow.
Lemma 3.1. Denote by Z the compact set of R n representing the spatial domain of the system. Then denote U, V two sets of smooth functions from Z to R q and R p respectively. Consider now a matrix differential operator L and denote by L * its formal adjoint ([7] Definition 3.2). Then, for every function u ∈ U and v ∈ V ( with u, v, Lu and L * v all in L 2 (Z)), there exists a differential operator B L defined on the boundary ∂Z (B L (u, v) ∈ L 2 (∂Z)) such that: We say that B L is the differential operator induced on the boundary ∂Z by the differential operator L.
For example, if we consider the first order spatial derivative L = ∂ z on the domain Z = [a, b], forall u, v ∈ Z, then we get: Therefore the boundary operator B ∂z induced by L = ∂ z is here simply the evaluation of the inner product u, v = uv on the boundary ∂Z = {a, b}. Such induced boundary operators may be constructed for a larger class higher-order skew symmetric differential operators. The resulting class of systems of the form (3.1) together with the class of input-output variables which are admissible in order to define a well-posed linear problem and generate a contraction semigroup are defined in [5,14]. It must be noticed that the class of systems of the form (3.1) includes most classical hyperbolic examples such as wave, membrane, plates and beams equations, shallow water, Boussinesq, Korteg de Vries and Navier-Stokes flow equations, Maxwell field equations, etc. but also some parabolic examples as it will be shown hereafter with the plasma poloidal flux resistive diffusion equation.
The purpose of the feedback control to be designed is to match a target (or "desired") canonical passive port Hamiltonian system of the form: . and H d are respectively the desired system interconnection operator, damping operator and hamiltonian. B d denotes a desired boundary operator which ensures that the couple (ỹ,ũ 2 ) is a passive input-output pair for the target system with respect to the storage functional H d . Again the definition and parametrization of all admissible passive input-output pairs which lead to a well-posed boundary control system may be found in [6] or [15,16]. Note that the passivity with the so-called impedance-passive input-output pairs of variables simply results from the energy balance equation which reads: (3.6) thanks to the skew-symmetry of the differential operator J d 1 . B d is the differential operator induced on the boundary ∂Z = {a, b} by the skew-adjoint operator J d (cf. [7] Corollary 3.2 or lemma 3.1 ).

Remark 1.
• Unlike the "traditional" approach which consists in homogenizing the boundary control to embed it in the system state equation and then handle it as a distributed control, the method presented in this paper reverses the idea. In other words, the distributed control is used to transform the original system (3.1) into the target canonical boundary control port-Hamiltonian system (3.5) with a boundary control (usually not the same as the one in the original system). • A dissipative target port-Hamiltonian system, that is with R d > 0, is asymptotically stable even with an homogenous boundary conditionũ 2 = 0. Otherwise, when R d ≥ 0 is only positive semidefinite, a simple "boundary damping" injection of the formũ 2 = −K pỹ2 , K p > 0 ensures the stabilization (cf. [14]). The matching equation then determines the "distributed" control gu 1 .

3.2.
Energy shaping and damping assignment for a subclass of linear port-Hamiltonian systems. We focus in this paper on a subclass of linear port-Hamiltonian systems. Assume that the energy function of the original system is the quadratic form H =´Z 1 2 x T QxdV (which means that this original system is linear) and that the canonical choice (u 2 , y 2 ) = (∂ x H) | ∂Z of boundary input-output impedance-passive variables has been selected [14]. It may happen that the interconnection operator J has a suitable form which should remain unchanged in the target system. In fact, this is the most common case and the geometric interconnection structure of the actual model should not be changed in the target system unless some specific purpose is given since it affects the structural invariants and intrinsic dynamical behaviour of the system. For instance, in the resistive diffusion equation example hereafter, the interconnection operator J is defined using the canonical derivation operator together with some given boundary conditions (see section 3) and should be preserved since it implies (with a dissipative closure equation) a purely dissipative input-output operator with a spectrum entirely lying on the negative real half-axis in the complex plane. In such cases where the interconnection structure of the actual system must remain unchanged, the target system is obtained by using only "energy shaping" and/or "damping injection": The new "passive" boundary control in the target system is determined via the new Hamiltonian H d =´Z H d dV and thus may be related with the original system boundary control u 2 as: The new boundary controlũ 2 in the target system is modified only by the energy shaping H a and not influenced by the damping injection R a . Furthermore, in the considered linear case, exponential stability is achieved even without any supplementary boundary controlũ 2 . The stability of the target system may be proved via the first and second Arnold's stability theorems [12] using some particular norm. In our case, choosing the norm associated with the energy stored in the target system results in very simple calculations to prove the asymptotic (exponential) stability with respect to this norm. Indeed, assume that the energy function of the target system is the quadratic form Considering now the case without boundary control (i.e.ũ 2 = 0), the energy balance (3.6) becomes: where Γ is the minimum singular value of (Q d ) In the case of a non-zero damping injection R d > 0 , we have Γ > 0 and: which proves the exponential stability in the sense of Lyapunov with respect to the energy norm . Q d . When R d = 0, exponential stability is not achieved but it is still possible to add a supplementary boundary damping of the formũ 2 = −K pỹ2 , K p > 0 (see [14]) to increase the convergence speed of the closed-loop system.

Average matching equation solutions.
How to solve the matching equation (3.7) and how to parametrize the solutions are major concerns in the IDA-PBC literature even for the control design for finite dimensional systems [9]. Since in our case the distributed control is only finite rank, there is no solution in the general case for the matching equation in the infinite dimensional case. We present hereafter two approaches to "solve" this problem.
• A first case occurs when there exists an adjoint function g ⊥ of g so-that their inner product g ⊥ g = 0, and once two among three control parameters J d , R d , and H d are set. Then the third target parameter is the solution of the linear equation: Note that the parameters defined in this method should respect the structural It is not always feasible to guarantee the existence of a solution with these properties for the equation here above. However, it exists for the resistive diffusion equation in the case of some damping assignment control as it will be shown here after in section 3. The idea is thus to restrict sufficiently the class of admissible target systems.
• A second approach consists in solving the matching equation only in an average sense. Indeed the scalar value of u 1 (t) not depending on the spatial coordinate z may be isolated from the control spatial distribution g (z, t) in the matching equation (3.13) and thus may be extracted as (3.14) • Therefore the obtained control u 1 is the one which cancels the average value of the residual for the matching equation on whole spatial domain [0, 1]. This idea could be extended to higher moments of the matching error residual in the case were several control variables are available.

Remarks.
• The previous solutions for the matching equation impose the errors from the average calculus. Therefore, one can't ensure the exact transformation from the original system into the target one, which is stable, with such derived control. The discussion then consists in knowing whether or not the boundary controlũ 2 can be used to stabilize this "matching error". • If one can prove that there exists a solution u 1 to the matching equation (3.7) which can bring the system to the stable target form (3.5), the convergence rate can be further improved using the boundary controlũ 2 . In other words, it is possible to determine a feedback control u 1 (t) which guarantees the existence of a solution to the matching equation without expliciting the choice of the control parameters J d , R d , and H d but only guarantees their existence. The example of a damping injection control design for the resistive diffusion equation will be cast using this approach in the next section.
The resistive diffusion of plasma poloidal flux in Tokamak will be figured out in the next part as an example of the previous proposition. Then the boundary control is used to accelerate the convergence. In the next subsection we recall and revisit the port Hamiltonian formulation for the resistive diffusion equation [17,18]. Then we will present the controller tuning method, some simulation results and robustness of the corresponding controller with respect to errors on some parameter (the resistivity) and on the distributed control (the non inductive current) in the following subsections.

4.1.
Infinite-dimensional PCH formulation for the resistive diffusion equation. The resistive diffusion equation for the poloidal magnetic flux [2] may be expressed in canonical port Hamiltonian form using some specific 1D variables var which depend only on the reduced radial coordinate ρ in a toroidal coordinate systems (see figure 4.1). In fact in these toric magnetic coordinates, the spatial coordinate z is in fact the normalized index of a corresponding nested magnetic surface (those surfaces all have basically the form of a deformed torus). The components of all vector fields in the radial ρ direction are zeros since these magnetic surfaces may be proven to be isobaric, isothermic and iso-flux at the same time [2]. The port-Hamiltonian formulation of the resistive diffusion equation in these  Magnetic toric coordinates (ρ, θ, φ). B θ and B φ are two magnetic field coordinates; B ρ = 0; R 0 denotes the tokamak major radius and I p is the total plasma current magnetic toric coordinates reads: In this model we will focus on the specific flow variables ∂ t (−zD φ ), ∂ (−R 0 B θ ) which are used to compute the poloidal magnetic flux and correspond respectively to the electric intensity and the magnetic density flows. The total current density J ext + J bs includes the bootstrap current J bs described in [21] (a magnetohydrodynamic coupling effect which produces and extra current density) and the external current source J ext whose the shape function is f ext (z) (depending on actuator properties) so-that J ext = f ext P ext . R 0 is the major plasma radius while z is the spatial index between 0 and 1 , from the center to the normalized minor plasma radius (corresponding to the last closed magnetic surface). C 2 (z) , C 3 (z) are the coordinate coefficients defined in [2]. The plasma resistivity η is a parameter varying significantly with the plasma temperature [22]. However in this paper we do not consider Thermo-Magneto-Hydro-Dynamical couplings. Hence the plasma resistivity is simply considered as a state and time dependent variable η (z, t). The electric permittivity and magnetic permeability µ are considered to be the void permittivity and void permeability since tokamaks are operating at very low densities. The boundary variables correspond to the total plasma current I p1 and the loop voltage V loop produced by the external electric coils. Notice also that there's no injection energy source at the center of the tokamak, i.e. f ∂0 e ∂0 = 0. The chosen Hamiltonian function for this model is the usual electromagnetic energy quadratic form: with Qx = ∂ x H as given in (4.1). Adopting notations from previous sections, this port-Hamiltonian model may be written in the form: Remark 3. In general, the control gu 1 is not "fully distributed", it is regulated only by the scalar power u 1 (t), g (x, t) is a function of system state and time. A feedforward computation is mandatory in order to determine the accessible steady state with respect to the actuator constraint g. Thus the error system (4.3) is in fact the one defined for the error after linearization around the chosen equilibrium.
Remark 4. The original system (4.3) without distributed control u 1 makes uses of the well known canonical Stokes-Dirac interconnection structure [13] and it is easily seen that its solution (if any exists) must satisfy the passivity property: The system dissipation R ≥ 0 is not strictly positive definite and the distributed control gu 1 (and the boundary control) will be used therefore to get asymptotic stability and improve the convergence speed.

4.2.
Controller tuning. We will first consider a damping assignment for the error system in (4.3), using only the distributed control to modify the system dissipation, while preserving the stored energy and Q d = Q (i.e. Q a = 0, no energy shaping). We thus want to perform the matching: Denote now R d = r 1 r 12 r 12 r 2 > 0, the desired strictly positive definite damping operator with: and define R 1 = C3 η ; Q = In this simplest case, the matching of the distributed controller term gu 1 which ensures the desired dissipation R d reads: Therefore, one only has to prove the point wise existence of this scalar control u 1 (t) for all x value. Meanwhile, the boundary control u 2 will be designed to accelerate the convergence of the solution. The matching equation (4.7) leads to: One can easily derive from (4.6) and (4.9) the solution: with the constraints or "bounds" on u 1 : The constraints in (4.11) does not guarantee the existence of a solution u 1 since the term −f ext Q 1 x 1 , depending on the system state x 1 , has undefined sign and value.
However, in the special case where r 1 > R 1 , u 1 = 0 is a trivial solution. There are then other values of u 1 satisfying the constraints. The original system then becomes the target system with the energy balance in (3.6). Moreover, a supplementary boundary damping injectionũ 2 = −K pỹ2 , K p > 0 will accelerate the convergence of the closed loop system to the equilibrium 0 (that is the convergence of the real system to the desired state used to design the feedforward control).

Simulation result.
The above IDA-PBC control law is applied to the resistive diffusion model in Tokamak plasma. However, this model doesn't stand alone, since it requires the time-variant profiles of resistivity and bootstrap current. The plasma resistivity mainly depends on the temperature whereas the bootstrap current is mainly function of its gradient. The Tore-Supra WEST Tokamak configuration is considered with the plasma parameters given as in [4].
Two separated step references of I p1 at t = 7s and I p0 (the plasma current at the center z ≈ 0) at t = 11s are used to illustrate the behavior of this control law. First of all, since the system is naturally dissipative, the feedforward in figure 4.2 shows the convergence of the open-loop system. Thanks to the re-computation of equilibrium at each step time, the desired reference is reached at the steady state. Then in figure 4.3 the response time is decreased by the feedback effect via u 1 . The convergence speed does improve as expected. In fact, the tuning of u 1 is not free since (4.11) must hold. The results with the complete feedback control including the finite rank distributed control u 1 and the boundary control u 2 are shown in the figure 4.4. The boundary plasma I p1 quickly reaches the reference as the boundary control V loop affects directly the conjugated variable I p1 . The boundary effect needs more time to propagate to the center, in order to converge I p0 .  Remark 5. The robustness analysis of the controller respects to two kinds of uncertainties: the major uncertainties (and model errors) are those on the system dissipation R resulting from poor estimations for the plasma resistivity η (z, t) and uncertainties related to the linearization assumption made in the derivation the feedforward control and in the estimation of the bootstrap current J bs (z, t) (the bootstrap current is taken account in the feedforward generation). According to [1], the system will remain stable as long as the desired dissipation is still positive and high enough to compensate the bounded perturbations. A "strong" enough  4.4. Robustness analysis. The robustness of the controller with respect to two kinds of uncertainties is investigated. In the resistive diffusion example, the major uncertainties (and model errors) are those on the system dissipation R resulting from poor estimations for the plasma resistivity η (z, t) and uncertainties related to the linearization assumption made in the derivation the feedforward control and in the estimation of the bootstrap current J bs (z, t) (the bootstrap current is taken account in the feedforward generation). Uncertainties on the resistivity. They lead to the parameter uncertainties δR = δR T applied on the dissipation matrix R. The disturbed system can then be defined as: 12) The perturbed system will remain stable as long as the total dissipation remains positive, i.e. [R d + δR] > 0. Uncertainties on the bootstrap estimation and from the linearization. These uncertainties are represented in the form of a disturbance ζ in the closed-loop system: ∂H d ∂x (4.14) (see [1] for explanations on this condition of disturbance rejection for IDA-PBC controllers) . This robustness analysis is important for the considered example since there are simulateneous uncertainties on the resistivity and on the bootstrap currents which both strongly depend on the temperature profiles. The figures 4.5 and 4.6 show the result when these disturbances are taken into account. The differences between the new actions and the one in the ideal case without disturbances are there denoted (δP ext , δV loop ). In figure 4.5 , it may be observed that the closed-loop system is still stable but with a static error. An integrator is then used to cancel this static error. The ideas behind the inclusion of an integrator effect in the IDA-PBC design are from [10] and are here tested on this infinite-dimensional example. In figure 4.6, it may be observed that the feedback control with the supplementary integrator compensates the perturbations and slowly recovers the correct values for the equilibrium. The plasma current at the center I p0 always requires more time to converge to the reference value (than I p1 , the one at the external plasma boundary) since the boundary control action is at the external boundary and the progressively diffuse troughout the spatial domain.  Simulation results with uncertainties on η (5% from t = 6s) and J bs (20% from t = 12s).

Conclusion.
A damping assignment control has been designed for a finite rank distributed control for the resistive diffusion model of the poloidal magnetic flux in tokamak reactors. The control design is based on the port Hamiltonian formulation of the resistive diffusion equation. The actuator spatial distribution has been included in the control design and the matching of the resulting controlled system with the asymptotically stable target system is guaranteed. This is not the case when we use some average matching method for the infinite dimensional system or any finite dimensional IDA-PBC controller.
The proposed controller has been tested on a simulation tool developed for the Tore Supra WEST configuration (experimental facility at CEA, Cadarache). Numerical experiments show that indeed asymptotic convergence is reached with this feedback control. Moreover the controller is shown to be robust against the errors on the resistivity and on the non inductive current radial shape which are the two major physical uncertainties in the resistive diffusion example.
In this preliminary work, the controller parameters choices for the target system have not been optimized and much work remains to determine how to choose J d , R d , and H d in order to minimize the residue in the matching equation in the general case and how to stabilize the matching error with the boundary controlũ 2 . Besides, damping assignment has been realized by using the distributed finite rank control u 1 only and modifying the interconnection structure as well as shaping the  energy of the system should still be realized with the controlũ 2 . before one can talk about a real IDA-PBC method.