Boundary energy-shaping control of an isothermal tubular reactor

ABSTRACT

This paper illustrates a general synthesis methodology of asymptotic stabilizing, energy-based, boundary control laws that are applicable to a large class of distributed port-Hamiltonian systems. The methodological results are applied on a simplified model of an isothermal tubular reactor. Due to the presence of diffusion and convection, such example, even if relatively easy from a computational point of view, is not trivial. The idea here is to design a state feedback law able to perform the energy-shaping task, i.e. able to render the closed-loop system a port-Hamiltonian system with the same structure, but characterized by a new Hamiltonian with a unique and isolated minimum at the equilibrium. Asymptotic stability is then obtained via damping injection on the boundary and is a consequence of the LaSalle’s Invariance Principle in infinite dimensions.

KEYWORDS:

• Distributed port-Hamiltonian systems
• boundary control
• energy-shaping control
• stability of PDEs
• irreversible thermodynamics

1. Introduction

Port-Hamiltonian systems have been introduced about 20 years ago as the mathematical formalization of bond-graphs to describe lumped parameter physical systems in a unified manner [1–4]. The generalization to the infinite dimensional scenario leads to the definition of distributed port-Hamiltonian systems (see e.g. [5,6]) that have been introduced about one decade ago, and that have proved to represent a powerful framework for modelling, simulation and control of physical systems described by PDEs. Distributed port-Hamiltonian systems share analogous geometric properties with their finite dimensional counterpart, and also the development of stabilizing control laws follows the same rationale of the lumped parameter case. Since in most of the cases the Hamiltonian is the total energy of the system, stabilization could be obtained by driving the Hamiltonian to zero. As a consequence, having such a physical quantity at our disposal simplifies the controller design considerably.

Most of the current research on the stability and stabilization of distributed port-Hamiltonian systems deals with the development of boundary controllers. For example, in [7–15], this task has been accomplished by looking at, or generating, a set of Casimir functions in closed-loop that robustly (i.e. independently from the Hamiltonian functions) relates the state of the infinite dimensional port-Hamiltonian system with the state of the controller, which is a finite dimensional port-Hamiltonian system interconnected to the boundary of the distributed parameter one. The shape of the closed-loop energy function is changed by acting on the Hamiltonian of the controller. As discussed in [2,16,17], this procedure is the generalization of the control by interconnection via Casimir generation (energy-Casimir method) developed for finite dimensional systems. The result is an energy-balancing passivity-based controller that is not able to deal with equilibria that require an infinite amount of supplied energy in steady state, i.e. with the so-called dissipation obstacle.

In this paper, it is shown how to enlarge the class of boundary energy-shaping controllers by determining a feedback law that renders the closed-loop system a port-Hamiltonian system with the same structure as the open-loop system (i.e. Dirac and resistive structures are not modified), but characterized by a shaped Hamiltonian with the desired performances. This is the same concept adopted in finite dimensions in case of stabilization with state-modulated sources discussed in [16], or with the more general IDA-PBC control technique presented in [18]. The focus here is on the class of dissipative boundary control systems presented in [19], which generalizes the class of linear, distributed, port-Hamiltonian systems presented in [20,21] to include distributed dissipation and/or diffusive effects. Our main contribution is the following: by transforming the original system via state feedback into a new one with a Hamiltonian function that has an isolated minimum at the equilibrium, simple stability is obtained. Then, asymptotic stability is achieved by introducing damping by means of an additional control loop. In this respect, another important contribution of this paper is to show that, if it is possible to impose full boundary dissipation to the port-Hamiltonian system resulting from the energy-shaping procedure, then the desired equilibrium can be proved to be asymptotically stable.

The general methodology is applied to the boundary stabilization of an isothermal tubular reactor, but the proposed approach can be naturally extended to deal with the linearized model of non-isothermal tubular reactors with convection, dispersion and reaction phenomena. In this paper, the considered reaction is of type $A\to B$, and the system is modelled by the mass balance equation of element A given as a second-order linear PDE. Preliminarily, existence of solutions for this system in case of constant parameters, and a parametrization of all the input defining a boundary control system are discussed. Then, asymptotic stabilization of non-zero equilibria via energy shaping and damping injection is achieved in case the boundary input is associated to Dankwert conditions [22]. The presence of diffusion and convection makes this system not trivial and, at the best of our knowledge, it is the first time that the stabilization of systems in which such phenomena are present has been carried out in the port-Hamiltonian framework using a combination of energy-shaping and damping injection.

The paper is organized as follows. In Section 2, it is shown how the differential operators can be parametrized, and how from this parametrization one can define a boundary control system. In Section 3, the model of the tubular reactor and its associated boundary control system are presented. The energy-shaping boundary control technique together with asymptotic stability in the case of full boundary dissipation (damping injection) is discussed in Section 4. Then, in Section 5, the general methodology is applied to the boundary stabilization of the 1D tubular reactor, while conclusions and a discussion about possible future research activities are reported in Section 6.

2. Boundary control systems

We consider the class of dissipative boundary control systems:

(1) $\left\{\begin{array}{cc}\frac{\mathrm{\partial }x}{\mathrm{\partial }t}\left(t,z\right)& =\left[\mathcal{J}-\mathcal{G}S\left(z\right){\mathcal{G}}^{\ast }\right]\mathcal{L}\left(z\right)x\left(t,z\right)\\ u\left(t\right)& =\mathcal{B}x\left(t,z\right)\\ y\left(t\right)& =\mathcal{C}x\left(t,z\right)\end{array}\right)x\left(0,z\right)=x_0\left(z\right)$(1)

where $z\in [a,b]$ is the spatial coordinate, $S(\cdot )$ and $\mathcal{L}(\cdot )$ are bounded self-adjoint coercive operators on $L^2\left(\left(a,b\right);{\mathbb{R}}^n\right)$ that can depend on z, and $\mathcal{J}$ and $\mathcal{G}$ are differential operators. The state space is defined as $X\in L^2\left(\left(a,b\right);{\mathbb{R}}^n\right)$ with inner product ${⟨x_1,x_2⟩}_{\mathcal{L}}={⟨x_1,\mathcal{L}{x}_2⟩}_{\mathcal{L}}$ and norm ${∥x∥}_{\mathcal{L}}^2={⟨x,x⟩}_{\mathcal{L}}$. Hence, X is a Hilbert space. Note with this definition, the natural norm on X and the $L^2$ norm are equivalent. The reason for selecting this space is that $\parallel x{\parallel }_{\mathcal{L}}^2$ is usually equal to the energy function of the system and from a physical point of view, x and $\mathcal{L}x$ denote, respectively, the energy and co-energy variables of the system. Indeed, the energy of the system is given by the Hamiltonian

(2) $H\left(x\left(t,z\right)\right)=\frac{1}{2}\underset{a}{\overset{b}{\int }}{x}^{\mathrm{T}}\left(t,z\right)\mathcal{L}\left(z\right)x\left(t,z\right)\mathrm{d}z.$(2)

The differential operator $\mathcal{J}$ is formally skew-adjoint and is related to energy conservation and structural invariants. It is assumed in the form

(3) $\mathcal{J}x=\sum _{i=0}^N{P}_i\frac{{\mathrm{\partial }}^{i}x}{\mathrm{\partial }{z}^i}$(3)

with each $P_i$ an n × n real and constant matrix such that $P_i=(-1{)}^{i+1}{P}_i^{\mathrm{T}}$, $i=1,\dots ,N$. The differential operator $\mathcal{G}$ and its (formal) adjoint ${\mathcal{G}}^{\ast }$ are related to energy dissipation and defined as

(4) $\mathcal{G}x=\sum _{i=1}^N{G}_i\frac{{\mathrm{\partial }}^{i}x}{\mathrm{\partial }{z}^i}{\mathcal{G}}^{\ast }x=\sum _{i=1}^N{\left(-1\right)}^i{G}_i^{\mathrm{T}}\frac{{\mathrm{\partial }}^{i}x}{\mathrm{\partial }{z}^i}$(4)

with $G_i$ an n × m real and constant matrix. The operators $\mathcal{B}$ and $\mathcal{C}$ are, respectively, the boundary input and output operators, defined later on. System (1) can be rewritten using an extended skew-symmetric operator ${\mathcal{J}}_e$ given hereafter in terms of conjugate port variables, efforts $\left(e_S,e_R\right)$ and flows $\left(f_S,f_R\right)$ of proper dimensions, as follows:

(5) $\left(\begin{array}{c}{f}_S\\ f_R\end{array}\right)=\underset{=:{\mathcal{J}}_e}{\underbrace{\left(\begin{array}{cc}\mathcal{J}& \mathcal{G}\\ -{\mathcal{G}}^{\ast }& 0\end{array}\right)}}\left(\begin{array}{c}{e}_S\\ e_R\end{array}\right)e_R=S{f}_R$(5)

with $f_S(t,z)=\left(\mathrm{\partial }x/\mathrm{\partial }t\right)(t,z)$, $e_S(t,z)=\mathcal{L}(z)x(t,z)$ and $\left(f_R,e_R\right)$ the conjugated port variables associated with dissipation. The operator ${\mathcal{J}}_e$ makes explicit the interconnection structure associated with the power continuous energy flows in the system among the energy storage elements interconnected at the $\left(f_S,e_S\right)$ port and the dissipative elements at the $\left(f_R,e_R\right)$ port [6]. Notice that, just as within the framework of bond graphs, the effort variables are related to co-energy variables, while the flow variables to the time derivative of the energy variables. Physically, they usually correspond to mechanical–thermodynamic driving forces (intensive variables) and flows (time changes) of extensive variables, respectively.

Proposition 2.1

The operator ${\mathcal{J}}_e$ defined in (5), with $\mathcal{J}$ and $\mathcal{G}$ given in (3) and (4), respectively, is formally skew-symmetric and can be written as

(6) ${\mathcal{J}}_e\left(\begin{array}{c}{e}_S\\ e_R\end{array}\right)=\sum _{i=0}^N\underset{=:{\tilde{P}}_i}{\underbrace{\left(\begin{array}{cc}{P}_i& G_i\\ {\left(-1\right)}^{i+1}{G}_i^{\mathrm{T}}& 0\end{array}\right)}}\frac{{\mathrm{\partial }}^i}{\mathrm{\partial }{z}^i}\left(\begin{array}{c}{e}_S\\ e_R\end{array}\right)$(6)

with ${\tilde{P}}_i=(-1{)}^{i+1}{\tilde{P}}_i^{\mathrm{T}}$. Note that ${\tilde{P}}_N$ can have a rank deficiency.

Proof. See [20].□

In the following, we define the boundary port variables associated with the extended operator ${\mathcal{J}}_e$. A detailed derivation of the boundary port-variables can be found in [20,21], and in this section, we limit ourselves to introduce the technical tools that are later going to be used in the control design. To this end, let us first define the following matrices:

(7) $\tilde{Q}=\left(\begin{array}{ccccc}{\tilde{P}}_1& {\tilde{P}}_2& \cdots & {\tilde{P}}_{N-1}& {\tilde{P}}_N\\ -{\tilde{P}}_2& -{\tilde{P}}_3& \cdots & -{\tilde{P}}_N& 0\\ {\tilde{P}}_3& {\tilde{P}}_4& \cdots & 0& 0\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ {\left(-1\right)}^{N-1}{\tilde{P}}_N& 0& \cdots & 0& 0\end{array}\right)$(7)

and M spanning the column of $\tilde{Q}$ (used to take into account the possible rank deficiency of $\tilde{Q}$), i.e. ${\tilde{Q}}_1=M^{\mathrm{T}}\tilde{Q}M$, and $M_Q={\left(M^{\mathrm{T}}M\right)}^{-1}{M}^{\mathrm{T}}$. Then

(8) $R_{\mathrm{e}\mathrm{x}\mathrm{t}}=\frac{1}{\sqrt{2}}\left(\begin{array}{cc}{\tilde{Q}}_1& -{\tilde{Q}}_1\\ I& I\end{array}\right)\left(\begin{array}{cc}{M}_Q& 0\\ 0& M_Q\end{array}\right)$(8)

where I is the $(n+m)N\times (n+m)N$ identity matrix. Notice that $P_0$ does not play a role in the energy exchange at the boundaries, and then in the definition of the associated boundary port. The previous matrices follow from the fact that any skew-symmetric differential operator gives rise to a symmetric bilinear form on the space of boundary variables, where the coefficients of the operator are captured in the matrix $\tilde{Q}$. Then, $R_{\mathrm{e}\mathrm{x}\mathrm{t}}$ follows from the factorization of $\tilde{Q}$ in such a way that it allows to define the port-variables and at the same time the bilinear form becomes independent of the coefficients of the operator ${\mathcal{J}}_e$. Let us now define the boundary port variables associated with the differential operator ${\mathcal{J}}_e$, which are the vectors $e_{e,\mathrm{\partial }},\hspace{0.17em}{f}_{e,\mathrm{\partial }}\in {\mathbb{R}}^{2nN}$ defined by

(9) $\left(\begin{array}{c}{f}_{e,\mathrm{\partial }}\\ e_{e,\mathrm{\partial }}\end{array}\right)=R_{\mathrm{e}\mathrm{x}\mathrm{t}}\left(\begin{array}{c}{e}_e\left(b\right)\\ ⋮\\ \frac{{\mathrm{\partial }}^{N-1}{e}_e}{\mathrm{\partial }{z}^{N-1}}\left(b\right)\\ e_e\left(a\right)\\ ⋮\\ \frac{{\mathrm{\partial }}^{N-1}{e}_e}{\mathrm{\partial }{z}^{N-1}}\left(a\right)\end{array}\right)$(9)

where

(10) $e_e\left(t,z\right)=\left(\begin{array}{c}{e}_S\left(t,z\right)\\ -S\left(z\right){\mathcal{G}}^{\ast }{e}_S\left(t,z\right)\end{array}\right)=\underset{=:{\mathcal{G}}_e}{\underbrace{\left(\begin{array}{c}I\\ -S\left(z\right){\mathcal{G}}^{\ast }\end{array}\right)}}{e}_S\left(t,z\right)$(10)

is the extended vector of the effort variables. The definition of the boundary port-variables depends entirely on the coefficients of the operator ${\mathcal{J}}_e$. These port-variables should be seen as an operator acting on the boundary of the spatial domain. After using the definition of the boundary port-variables of $R_{\mathrm{e}\mathrm{x}\mathrm{t}}$ and of $\tilde{Q}$, it is possible to relate the integral over an interval to the values at the boundaries. Indeed, by using integration by part, we have

$\underset{a}{\overset{b}{\int }}\left(e_S^{\mathrm{T}}{f}_S+e_R^{\mathrm{T}}{f}_R\right)\mathrm{d}z=f_{e_S,\mathrm{\partial }}^{\mathrm{T}}{e}_{e_R,\mathrm{\partial }}+e_{e_S,\mathrm{\partial }}^{\mathrm{T}}{f}_{e_R,\mathrm{\partial }}$

The next theorem provides the matrix condition that has to be satisfied to ensure the existence of solutions for the PDE (1) and to define a boundary control system.

Theorem 2.2

(Boundary control systems): Let W be an $(n+m)N\times 2(n+m)N$ matrix. If W has full rank and satisfies the inequality:

(11) $W\mathrm{\Sigma }{W}^{\mathrm{T}}\ge 0$(11)

where

$\mathrm{\Sigma }=\left(\begin{array}{cc}0& I\\ I& 0\end{array}\right)$

then system (1) in which

(12) $u\left(t\right)=W\left(\begin{array}{c}{f}_{e,\mathrm{\partial }}\left(t\right)\\ e_{e,\mathrm{\partial }}\left(t\right)\end{array}\right)\equiv \mathcal{B}x\left(t,z\right)$(12)

is a boundary control system, and the operator $\mathcal{A}=\left(\mathcal{J}-\mathcal{G}S{\mathcal{G}}^{\ast }\right)\mathcal{L}$ with domain

(13) $D\left(\mathcal{A}\right)=\left\{e_S\in H^N\left(\left(a,b\right);{\mathbb{R}}^n\right)\left|S{\mathcal{G}}^{\ast }{e}_S\in H^N\left(\left(a,b\right);{\mathbb{R}}^n\right),\left(\begin{array}{c}{f}_{e,\mathrm{\partial }}\\ e_{e,\mathrm{\partial }}\end{array}\right)\in KerW\right)\right\}$(13)

where $H^N\left((a,b);{\mathbb{R}}^n\right)$ denotes the Sobolev space of order N, generates a contraction semigroup on X. Moreover, let $\tilde{W}$ be a full rank $(n+m)N\times 2(n+m)N$ matrix such that $\left(\begin{array}{cc}{W}^{\mathrm{T}}& {\tilde{W}}^{\mathrm{T}}\end{array}\right)$ is invertible and let P be given by

$P={\left(\begin{array}{cc}W\mathrm{\Sigma }{W}^{\mathrm{T}}& W\mathrm{\Sigma }{\tilde{W}}^{\mathrm{T}}\\ \tilde{W}\mathrm{\Sigma }{W}^{\mathrm{T}}& \tilde{W}\mathrm{\Sigma }{\tilde{W}}^{\mathrm{T}}\end{array}\right)}^{-1}$

If the output is defined as

(14) $y(t)=\tilde{W}\left(\begin{array}{c}{f}_{\mathrm{\partial }}\left(t\right)\\ e_{\mathrm{\partial }}\left(t\right)\end{array}\right)=:\mathcal{C}x$(14)

then, for $u\in C^2\left((0,\mathrm{\infty });{\mathbb{R}}^{n+m}\right)$, the following energy balance equation is satisfied:

(15) $\frac{1}{2}\frac{\mathrm{d}H}{\mathrm{d}t}\left(x\left(t,\cdot \right)\right)\le \frac{1}{2}{\left(\begin{array}{c}u\left(t\right)\\ y\left(t\right)\end{array}\right)}^{\mathrm{T}}P\left(\begin{array}{c}u\left(t\right)\\ y\left(t\right)\end{array}\right).$(15)

Proof. See [19].□

It should be remarked that the boundary port-variables (10) are the power conjugated variables associated to the energy balance. The boundary input (12) corresponds to the imposed boundary conditions, i.e. the physical efforts and flows imposed to the system at the boundary of the spatial domain. Hence, the matrix W defines the possible linear combinations of the boundary variables which guarantee that the solutions of boundary control system exist, and for physical systems, these linear combinations are related to physically admissible boundary conditions.

3. Boundary control system associated with the tubular reactor

We consider now a 1D linear tubular reactor of length L, in which convection, (axial) dispersion and reaction phenomena occur. The mass balance equation on specie A leads to the following PDE:

(16) $\frac{\mathrm{\partial }{\rho }_A}{\mathrm{\partial }t}\left(t,z\right)=D\frac{{\mathrm{\partial }}^2{\rho }_A}{\mathrm{\partial }{z}^2}\left(t,z\right)-v\frac{\mathrm{\partial }{\rho }_A}{\mathrm{\partial }z}\left(t,z\right)-k{\rho }_A\left(t,z\right)z\in \left[0,L\right]$(16)

where ${\rho }_A$ is the mass concentration of element A, $D>0$ is the dispersion coefficient, $v>0$ is the velocity and $k>0$ the reaction kinetics. Hereafter, a possible way of writing (16) in the general form (1) is given. This formulation differs and improves what has been already presented e.g. in [19]. The proposed formulation is canonical since the differential operators $\mathcal{J}$ and $\mathcal{G}$ do not depend on the system parameters. Moreover, the conditions of Theorem 2.2 are satisfied provided that D, v and k are positive.

The first step consists in determining the differential operator ${\mathcal{J}}_e$ and the dissipative relation in (5). So, let us assume that $N=1$, $n=1$ and $m=2$, and consequently $e_S\in L^2((0,L);\mathbb{R})$ and $e_R=\left(e_{R1},e_{R2}\right)\in L^2\left((0,L);{\mathbb{R}}^2\right)$. Hence, $\mathcal{J}=-\left(\mathrm{\partial }/\mathrm{\partial }z\right)$ and $\mathcal{G}=\left(\begin{array}{cc}\left(\mathrm{\partial }/\mathrm{\partial }z\right)& 1\end{array}\right)$, with Hamiltonian function

(17) $H\left({\rho }_A\left(t,z\right)\right)=\frac{1}{2}\underset{0}{\overset{L}{\int }}v{\rho }_A^2\left(t,z\right)\mathrm{d}z$(17)

which is related to the mass balance. The extended differential operator ${\mathcal{J}}_e$ in (5) is then given by

(18) $\left(\begin{array}{c}{f}_S\\ f_{R1}\\ f_{R2}\end{array}\right)=\left(\begin{array}{ccc}-\frac{\mathrm{\partial }}{\mathrm{\partial }z}& \frac{\mathrm{\partial }}{\mathrm{\partial }z}& 1\\ \frac{\mathrm{\partial }}{\mathrm{\partial }z}& 0& 0\\ -1& 0& 0\end{array}\right)\left(\begin{array}{c}{e}_S\\ e_{R1}\\ e_{R2}\end{array}\right)$(18)

while the dissipative relation, defined by S, by

(19) $\left(\begin{array}{c}{e}_{R1}\\ e_{R2}\end{array}\right)=\left(\begin{array}{cc}\frac{D}{v}& 0\\ 0& \frac{k}{v}\end{array}\right)\left(\begin{array}{c}{f}_{R1}\\ f_{R2}\end{array}\right).$(19)

We see that the conserved quantity, characterized by $\mathcal{J}$, is related to the convection phenomenon, while the dissipation, characterized by $\mathcal{G}$, is related to the diffusion and reaction phenomena. This is indeed coherent with previous existing energy-based models of tubular reactors [19]. It follows that

(20) $\tilde{Q}={\tilde{P}}_1=\left(\begin{array}{ccc}-1& 1& 0\\ 1& 0& 0\\ 0& 0& 0\end{array}\right)$(20)

which is not full rank. This implies that $\left(f_{R2},e_{R2}\right)$ does not play any active role in terms of the power exchange at the boundary for (16). One can define M as:

$M=\left(\begin{array}{cc}1& 0\\ 0& 1\\ 0& 0\end{array}\right)$

and the boundary port variables (9) are obtained with $R_{\mathrm{e}\mathrm{x}\mathrm{t}}$ in (8) in which

${\tilde{Q}}_1=M^{\mathrm{T}}\tilde{Q}M=\left(\begin{array}{cc}-1& 1\\ 1& 0\end{array}\right).$

This corresponds to a projection that allows to preserve the Dirac structure in case ${\tilde{P}}_N$ in (6), and then in (7), is not full rank. In this case, $N=1$. More details e.g. in [23]. This implies that

(21) $R_{\mathrm{e}\mathrm{x}\mathrm{t}}=\frac{1}{\sqrt{2}}\left(\begin{array}{cccc}-1& 1& 1& -1\\ 1& 0& -1& 0\\ 1& 0& 1& 0\\ 0& 1& 0& 1\end{array}\right)$(21)

and that

(22) $e_e=\left(\begin{array}{c}{e}_S\\ e_{R1}\end{array}\right)=\left(\begin{array}{c}{e}_S\\ \frac{D}{v}\frac{\mathrm{\partial }{e}_S}{\mathrm{\partial }z}\end{array}\right)=\left(\begin{array}{c}1\\ \frac{D}{v}\frac{\mathrm{\partial }}{\mathrm{\partial }z}\end{array}\right)e_S.$(22)

The PDE (16) follows from (18) and (19) provided that $f_S=\left(\mathrm{\partial }x/\mathrm{\partial }t\right)$ and $e_S=v{\rho }_A$, which in turns implies that $\mathcal{L}=v$. In terms of the state variable ${\rho }_A$, the boundary port variables are then given by

(23) $\left(\begin{array}{c}{f}_{e,\mathrm{\partial }}\\ e_{e,\mathrm{\partial }}\end{array}\right)=\frac{1}{\sqrt{2}}\left(\begin{array}{c}v\left({\rho }_A\left(0\right)-{\rho }_A\left(L\right)\right)-D\left(\frac{\mathrm{\partial }{\rho }_A}{\mathrm{\partial }z}\left(0\right)-\frac{\mathrm{\partial }{\rho }_A}{\mathrm{\partial }z}\left(L\right)\right)\\ -v\left({\rho }_A\left(0\right)-{\rho }_A\left(L\right)\right)\\ v\left({\rho }_A\left(0\right)+{\rho }_A\left(L\right)\right)\\ D\left(\frac{\mathrm{\partial }{\rho }_A}{\mathrm{\partial }z}\left(0\right)+\frac{\mathrm{\partial }{\rho }_A}{\mathrm{\partial }z}\left(L\right)\right)\end{array}\right).$(23)

Let us now consider the Danckwert boundary conditions [22] which consist to fix the total inlet flow $v{\rho }_A-D\left(\mathrm{\partial }{\rho }_A/\mathrm{\partial }z\right)$ at $z=0$, and to assume a vanishing outlet dispersion flow at $z=L$. This is achieved by choosing:

(24) $u=\left(\begin{array}{c}v{\rho }_A\left(0\right)-D\frac{\mathrm{\partial }{\rho }_A}{\mathrm{\partial }z}\left(0\right)\\ D\frac{\mathrm{\partial }{\rho }_A}{\mathrm{\partial }z}\left(L\right)\end{array}\right)$(24)

which can be written in the form (12) if

(25) $W=\frac{\sqrt{2}}{2}\left(\begin{array}{cccc}1& 0& 1& -1\\ 1& 1& 0& 1\end{array}\right).$(25)

With Theorem 2.2 in mind, this system defines a boundary control system, with contraction semigroup generator if condition (11) holds. With simple computations, we get that $W\mathrm{\Sigma }{W}^{\mathrm{T}}=I$, so (16) with boundary input given in (24) defines a boundary control systems. Note that, differently from [19] since the operators $\mathcal{J}$ and $\mathcal{G}$ are canonical, this property holds also when that parameters v, D and k depend on the z coordinate.

4. Boundary control by energy shaping

In this section, it is shown how to design a boundary state-feedback in the form

(26) $u\left(t\right)=\beta \left(x\left(t,\cdot \right)\right)+u^{\prime }\left(t\right)$(26)

which renders (1) into the closed-loop form

(27) $\left\{\begin{array}{cc}\frac{\mathrm{\partial }x}{\mathrm{\partial }t}\left(t,z\right)& =\left[\mathcal{J}-\mathcal{G}S\left(z\right){\mathcal{G}}^{\ast }\right]\frac{\delta H_d}{\delta x}\left(x\left(t,z\right)\right)\\ u^{\prime }\left(t\right)& ={\mathcal{B}}^{\prime }x\left(t,z\right)\end{array}\right)$(27)

where $H_d(x)=H(x)+H_a(x)$ is a desired closed-loop Hamiltonian function, that is properly designed in order to have a (possibly) global minimum at the equilibrium, and $\left(\delta /\delta x\right)$ denotes the functional (Frechét) derivative [5,24]. Note that, from (2), $\delta H/\delta x=\mathcal{L}x$. The target system (27) has the same internal structure of the original one, i.e. the operators $\mathcal{J}$ and $\mathcal{G}$, and the dissipative relation S are not changed, but a different Hamiltonian $H_d$, and boundary input $u^{\prime }$ are considered. For simplicity, we assume that there are no changes in the input–output port of the system, i.e. that $u^{\prime }$ is given by (9) and (12), but with respect to the desired closed-loop Hamiltonian, i.e.

(28) $e_e\left(t,z\right)={\mathcal{G}}_e\frac{\delta H_d}{\delta x}\left(x\left(t,z\right)\right)$(28)

where ${\mathcal{G}}_e$ is defined in (10), and W given as in (25).

The idea here is to overcome the intrinsic limitations of the energy-Casimir method that are associated with the admissible (internal) dissipation and to realize an ‘explicit’ energy-shaping procedure, as for the lumped parameter case. The energy-Casimir method, in fact, is constructive and it is based on the definition of a boundary controller in port-Hamiltonian form whose structure is chosen in order to have a proper set of Casimir functions in closed-loop. Such invariants show how it is possible to shape the total energy of the system by acting on the controller Hamiltonian: the stabilizing control law depends then on this choice and on the relation between states of the plant and of the controller specified by the Casimir functions. The associated control action can be written in terms of a state-feedback law $\beta \left(x\left(t,\cdot \right)\right)$ as in (27), but unfortunately, such control action suffers of all the intrinsic limitations of the energy-Casimir method, i.e. it is not able to stabilize equilibria that require an infinite amount of supplied energy in steady state [2,16]. This is the typical situation in systems with convective and dissipative effects, like (16).

Proposition 4.1

(Energy-shaping): Consider the boundary control system of Theorem 2.2 with Hamiltonian H defined in (2). Then, the boundary state-feedback law $u=\beta (x)+u^{\prime }$, being $u^{\prime }$ an auxiliary boundary input, renders (1) as the target dynamical system (27), with $H_d(x)=H(x)+H_a(x)$, if

(29) $\left[\mathcal{J}-\mathcal{G}S\left(z\right){\mathcal{G}}^{\ast }\right]\frac{\delta H_a}{\delta x}\left(x\right)=0,$(29)

(30) $\beta \left(x\right)+\left(\mathcal{B}-{\mathcal{B}}^{\prime }\right)x=0.$(30)

Proof. The proof is immediate by comparison of initial and target dynamics.□

With $u^{\prime }=0$, energy is not increasing along the trajectories of (27), i.e. ${\dot{H}}_d(x(t,\cdot ))\le 0$. Asymptotic stability can be then obtained by damping injection, provided that an output $y^{\prime }$ is properly defined. In this respect, with Theorem 2.2 in mind, the common choice (see e.g. [25–28].) is to have open-loop system (1) and target dynamic (27) in impedance form, i.e. to have the input and output pair (12) and (14) for which W and $\tilde{W}$ are such that

$W\mathrm{\Sigma }{W}^{\mathrm{T}}=\tilde{W}\mathrm{\Sigma }{\tilde{W}}^{\mathrm{T}}=0,W\mathrm{\Sigma }{\tilde{W}}^{\mathrm{T}}=I.$

With this choice, from (15), we obtain that ${\dot{H}}_d\le {y^{\prime }}^{\mathrm{T}}{u}^{\prime }$. On the other hand, the choice (24) for the boundary input u associated to (16), and then also for $u^{\prime }$ associated to (27), satisfies $W\mathrm{\Sigma }{W}^{\mathrm{T}}=I$. This suggests to define $y^{\prime }$ by selecting a matrix $\tilde{W}$ such that $\tilde{W}\mathrm{\Sigma }{\tilde{W}}^{\mathrm{T}}=-I$. As a consequence, for the target dynamic we have that

$\frac{\mathrm{d}}{\mathrm{d}t}{H}_d\left(x\left(t\right)\right)\le \frac{1}{2}{u^{\prime }}^{\mathrm{T}}\left(t\right)u^{\prime }\left(t\right)-\frac{1}{2}{y^{\prime }}^{\mathrm{T}}\left(t\right)y^{\prime }\left(t\right)$

and the associated system is said to be in scattering form. Such new boundary port $\left(u^{\prime },y^{\prime }\right)$ has now to be terminated over a dissipative element to obtain asymptotic stability of equilibria or just to improve the convergence rate:

(31) $u^{\prime }\left(t\right)=\mathrm{\Xi }{y}^{\prime }\left(t\right),I-{\mathrm{\Xi }}^{\mathrm{T}}\mathrm{\Xi }>0$(31)

with $\mathrm{\Xi }$ a $N(n+m)\times N(n+m)$ constant matrix. In what follows, we shall show that the energy-shaping control law defined in (29) and (30), combined with the damping injection relation (31) asymptotically stabilizes (1) at the equilibrium $(\mathcal{L}x{)}_{\star }\in H^N\left((a,b);\hspace{0.17em}{\mathbb{R}}^n\right)$, which is solution of

(32) $\left[\mathcal{J}-\mathcal{G}S\left(z\right){\mathcal{G}}^{\ast }\right]{\left(\mathcal{L}x\right)}_{\star }\left(z\right)=0.$(32)

The result is an application of the LaSalle’s Invariance Principle in infinite dimensions, see e.g. [29]. The first step is to determine how to choose $H_a$ so that (29) holds.

Proposition 4.2

The functions $H_a$ solutions of (29) are in the form $H_a(x)={\hat{H}}_a(\xi (x))$, with

(33) $\xi \left(x\left(t,\cdot \right)\right)=\underset{a}{\overset{b}{\int }}{\hat{\mathrm{\Phi }}}^{\mathrm{T}}\left(z\right)x\left(t,z\right)\mathrm{d}z$(33)

where $\hat{\mathrm{\Phi }}(z)=\left({\mathrm{\Phi }}_1(z),\dots ,{\mathrm{\Phi }}_{n_{\xi }}(z)\right)$. Here, the functions ${\mathrm{\Phi }}_i\in H^N\left((a,b);{\mathbb{R}}^n\right)$, $i=1,\dots ,n_{\xi }\le n$, are independent solutions of

(34) $\left[\mathcal{J}-\mathcal{G}S\left(z\right){\mathcal{G}}^{\ast }\right]{\mathrm{\Phi }}_i\left(z\right)=0.$(34)

Since in this paper we have restricted ourselves to the linear case, let us assume that $H_a$ is quadratic in $\xi$. Furthermore, denote by ${\varphi }_{\star }\in {\mathbb{R}}^{n_{\xi }}$ a vector such that $(\mathcal{L}x{)}_{\star }(z)=\hat{\mathrm{\Phi }}(z){\varphi }_{\star }$. Then, $H_d=H+H_a$ has a global minimum in $(\mathcal{L}x{)}_{\star }$ if

(35) $H_a\left(x\right)=\frac{1}{2}{\left\{\underset{a}{\overset{b}{\int }}{\hat{\mathrm{\Phi }}}^{\mathrm{T}}\left[x-{\mathcal{L}}^{-1}\hat{\mathrm{\Phi }}{\varphi }_{\star }\right]\mathrm{d}z\right\}}^{\mathrm{T}}\times Q_a\left\{\underset{a}{\overset{b}{\int }}{\hat{\mathrm{\Phi }}}^{\mathrm{T}}\left[x-{\mathcal{L}}^{-1}\hat{\mathrm{\Phi }}{\varphi }_{\star }\right]\mathrm{d}z\right\}-{\varphi }_{\star }^{\mathrm{T}}\left(\underset{a}{\overset{b}{\int }}{\hat{\mathrm{\Phi }}}^{\mathrm{T}}x\mathrm{d}z\right)+\kappa$(35)

where $Q_a={Q_a}^{\mathrm{T}}\ge 0$ and $\kappa \in \mathbb{R}$ is some constant.

Theorem 4.3

(Asymptotic stability): Let us consider the linear, infinite dimensional, port-Hamiltonian system (1) and the equilibrium $(\mathcal{L}x{)}_{\star }$ satisfying (32). Then, the control action $u=\beta (x)+u^{\prime }$ with $\beta$ defined in (30), being $H_a$ chosen as in (35), and with $u^{\prime }$ defined in (31), makes $(\mathcal{L}x{)}_{\star }$ asymptotically stable.

Proof. The closed-loop system subject to conditions (30)–(32) and (35) generates a contraction semigroup by Theorem 2.2. Thus, for any $x(0)\in X$, the solution $x(t)=T(t)x(0)$ (classical or weak) is bounded in X. Considering appropriate inner product and port variables the results of [30] (pages 49–50, Section 2.5.2) can be extended to operators ${\tilde{P}}_N$ with rank deficiency since, as discussed in [30], boundary inputs can be selected in such a way that the operator ${\mathcal{J}}_e$ defined in (6) has a compact resolvent and generates a contraction semigroup, similarly to what happens in the full rank case, treated in [23]. Hence $\left(I-{\mathcal{J}}_e\right)$, where I is the identity operator in X, which is a bounded invertible operator and from Theorem 2.26 ([30], page 49, Section 2.5.2), the resolvent ${\left(\lambda I-{\mathcal{J}}_e\right)}^{-1}$ is a compact operator for $\lambda >0$. From this, it follows that the trajectory of the solution $x(t)$, i.e. the set $\gamma (x(0))=\{x(t)\in X,t\ge 0\}$ is pre-compact in X (see Theorem 3.65 of [29]). The pre-compactness of the orbits of the closed-loop system implies asymptotic stability in spite of LaSalle’s Invariance Principle. In fact, let us assume for simplicity and without loss of generality that $\hat{\mathrm{\Phi }}(z)=(\mathcal{L}x{)}_{\star }(z)$, so that $\xi \in \mathbb{R}$ and $Q_a\in \mathbb{R}$. Select $Q_a>0$ and note that with a simple change of coordinates studying the stability of $(\mathcal{L}x{)}_{\star }$ is equivalent to studying the stability of the origin. For the closed-loop system, the following energy-balancing relation holds true:

$\frac{\mathrm{d}}{\mathrm{d}t}{H}_d\left(x\right)=-\underset{a}{\overset{b}{\int }}{\left(\frac{\delta H_d}{\delta x}\right)}^{\mathrm{T}}\mathcal{G}S\left(z\right){\mathcal{G}}^{\ast }\frac{\delta H_d}{\delta x}\mathrm{d}z-{y^{\prime }}^{\mathrm{T}}\left(I-{\mathrm{\Xi }}^{\mathrm{T}}\mathrm{\Xi }\right)y^{\prime }\le 0$

with $y^{\prime }$ defined such that it only vanishes at the desired equilibrium. Since $I-{\mathrm{\Xi }}^{\mathrm{T}}\mathrm{\Xi }$ is non-singular, it is not difficult to verify that energy is decreasing until a steady-state configuration $\bar{x}(t,z)$ is reached. Such configuration, possibly time-variant, satisfies

(36) $\begin{array}{cc}\mathcal{G}S\left(z\right){\mathcal{G}}^{\ast }\left[\left(\mathcal{L}\bar{x}\right)\left(t,z\right)+Q_a\bar{\xi }\left(t\right)\hat{\mathrm{\Phi }}\left(z\right)\right]& =0\\ \left(\begin{array}{c}\frac{\delta H_d}{\delta x}\left(\bar{x}\left(t,b\right)\right)\\ \frac{\delta H_d}{\delta x}\left(\bar{x}\left(t,a\right)\right)\end{array}\right)=\left(\begin{array}{c}\left(\mathcal{L}\bar{x}\right)\left(t,b\right)\\ \left(\mathcal{L}\bar{x}\right)\left(t,a\right)\end{array}\right)+Q_a\bar{\xi }\left(t\right)\left(\begin{array}{c}\hat{\mathrm{\Phi }}\left(b\right)\\ \hat{\mathrm{\Phi }}\left(a\right)\end{array}\right)& =0\end{array}$(36)

where $\bar{\xi }$ is the corresponding steady-state evolution of $\xi$. The second relation in (36) is a consequence of Theorem 2.2. To prove asymptotic stability by relying on LaSalle’s Invariance Principle [29], it is necessary to verify that the only steady state solution $\bar{x}(t,z)$ which is invariant and compatible with ${\dot{H}}_d=0$ is the origin. In this respect, we see that $\dot{\overline{\xi }}(t)=0$ means that $\bar{\xi }(t)={\xi }_{\star }$ in steady state. With some further computations, it is proved that $\varphi (t,z){:}_{=}\bar{x}(t,z)+{\xi }_{\star }{Q}_a{\mathcal{L}}^{-1}(z)\hat{\mathrm{\Phi }}(z)=0$ for $t\ge \tau$, being $\tau$ sufficiently large, which implies that $\bar{x}(t,z)=-{\xi }_{\star }{Q}_a{\mathcal{L}}^{-1}(z)\hat{\mathrm{\Phi }}(z)$ when $t\ge \tau$. From (14), we have that

${\xi }_{\star }=\underset{a}{\overset{b}{\int }}{\hat{\mathrm{\Phi }}}^{\mathrm{T}}\left(z\right)\bar{x}\left(t,z\right)\mathrm{d}z={\xi }_{\star }{Q}_a\underset{a}{\overset{b}{\int }}{\hat{\mathrm{\Phi }}}^{\mathrm{T}}\left(z\right){\mathcal{L}}^{-1}\left(z\right)\hat{\mathrm{\Phi }}\left(z\right)\mathrm{d}z$

Since the integral term is greater than 0 because $\mathcal{L}(\cdot )>0$, we have that ${\xi }_{\star }=0$, and then that $\bar{x}(t,z)=0$ for $t\ge \tau$. Then, the zero solution is the only invariant solution compatible with ${\dot{H}}_d=0$, which turns out to be asymptotically stable based on LaSalle’s Invariance Principle considerations.□

5. Asymptotic stabilization of the 1D linear tubular reactor

Let us consider the 1D tubular reactor model (16) equipped with the boundary input u defined in (23). The scope of this section is to show how to design a boundary state-feedback control law able to asymptotically stabilize an equilibrium solution ${\rho }_{A\star }\in L^2((0,L);\mathbb{R})$ of (16). The first step is to determine $\beta$ in (26) that allows to shape the open loop Hamiltonian (17), which we recall here: $H\left({\rho }_A(t,z)\right)=\left(1/2\right){\int }_0^{L}v{\rho }_A^2(t,z)\mathrm{d}z$. The open-loop Hamiltonian serves as starting point for the stabilization. Indeed, since the control objective is to introduce a global minimum at ${\rho }_{A\star }$ it is enough to shape (17) such that in closed-loop, it has a global minimum at the desired equilibrium. The solution of (32) becomes for this case

(37) $D\frac{{\mathrm{d}}^2{\rho }_{A\star }}{\mathrm{d}{z}^2}\left(z\right)-v\frac{\mathrm{d}{\rho }_{A\star }}{\mathrm{d}z}\left(z\right)-k{\rho }_{A\star }\left(z\right)=0$(37)

such that for some given $u_{\star }\in {\mathbb{R}}^2$

(38) $u_{\star }=\left(\begin{array}{c}v{\rho }_{A\star }(0)-D\frac{\mathrm{d}{\rho }_{A\star }}{\mathrm{d}z}(0)\\ D\frac{\mathrm{d}{\rho }_{A\star }}{\mathrm{d}z}(L)\end{array}\right)$(38)

where the causality associated with the choice (23) on the boundary input has been taken into account.

By following the energy-shaping methodology presented in Section 4, the PDE (29) provides the admissible functions $H_a$, while (30) the associated boundary control actions. With Proposition 4.2 in mind, the admissible functions $H_a$ take the form $H_a\left({\rho }_A\right)={\hat{H}}_a\left(\xi \left({\rho }_A\right)\right)$, with $\xi ({\rho }_A(t,\cdot ))={\int }_0^L{\rho }_{A\star }(t,z){\rho }_A(t,z)\hspace{0.17em}\mathrm{d}z$. A possible choice for $H_a$ is with

(39) ${\hat{H}}_a\left(\xi \right)=\frac{1}{2}K{\left(\xi -{\xi }_{\star }\right)}^2-v\left(\xi -{\xi }_{\star }\right)$(39)

where K is a positive gain and ${\xi }_k={\int }_0^L{\rho }_{A\ast }^2(z)dz$ the value of $\xi$ at the equilibrium. With this choice, the state feedback action $\beta$ obtained with (30) is able to shape the closed-loop Hamiltonian and to introduce a minimum in the desired equilibrium. Note that, when ${\rho }_A={\rho }_{A\star }$, we have $\beta \left({\rho }_a\right)=u_{\star }$ by construction. According to Theorem 4.3, asymptotic stability is obtained via damping injection (31) on the new control port $\left(u^{\prime },y^{\prime }\right)$, where $u^{\prime }$ is defined as

$u^{\prime }=\left(\begin{array}{c}\frac{\delta H_d}{\delta {\rho }_A}\left(0\right)-\frac{D}{v}\frac{\mathrm{\partial }}{\mathrm{\partial }z}\frac{\delta H_d}{\delta {\rho }_A}\left(0\right)\\ \frac{D}{v}\frac{\mathrm{\partial }}{\mathrm{\partial }z}\frac{\delta H_d}{\delta {\rho }_A}\left(L\right)\end{array}\right)$

with $H_d=H+H_a$, being H and $H_a$ as in (17) and (39), and

$y^{\prime }\left(t\right)=\left(\begin{array}{c}-\frac{\delta H_d}{\delta {\rho }_A}\left(L\right)+\frac{D}{v}\frac{\mathrm{\partial }}{\mathrm{\partial }z}\frac{\delta H_d}{\delta {\rho }_A}\left(L\right)\\ \frac{D}{v}\frac{\mathrm{\partial }}{\mathrm{\partial }z}\frac{\delta H_d}{\delta {\rho }_A}\left(0\right)\end{array}\right)$

which follows from (14) if $\tilde{W}=\left(\sqrt{2}/2\right)\left(\begin{array}{cccc}1& 0& -1& 1\\ -1& -1& 0& 1\end{array}\right)$.

6. Conclusions and future work

The motivating idea of the paper has been the development of a general synthesis methodology of control laws for a class of dissipative boundary control systems with 1D spatial domain. The feedback law is determined in such a way that its effect on the system is to shape the energy function, and to modify the dissipative structure at the boundary (damping injection). The first step is responsible for achieving simple stability of an equilibrium, while the second one for assuring asymptotic convergence of the trajectories. The problem of determining a feedback law able to shape the Hamiltonian in a proper manner has been tackled here by directly focusing on the trajectories of the open and closed-loop system, i.e. by determining the control action that renders the open-loop system into a new one, with the same structure but a different Hamiltonian in closed-loop. To achieve asymptotic stability, an additional feedback loop is implemented to introduce damping: the first loop obtained by applying the energy-shaping procedure is responsible for having a new Hamiltonian with an isolated minimum at the equilibrium, while the second one for dissipating energy until such minimum is reached. The resulting control law is proved to asymptotically stabilize the system. Such technique is applied to achieve asymptotic stability for a model of an isothermal tubular reactor with 1D spatial domain.

It is worth noting that, even if the proposed methodology has been developed for linear systems, some of the techniques discussed here can be generalized to cope with the non-linear case in a quite straightforward manner. This extension is the main future research topic, together with the stabilization of distributed port-Hamiltonian systems with 2D or 3D spatial domain.

Disclosure statement

No potential conflict of interest was reported by the authors.