Energy-optimal control of adaptive structures

Linear port-Hamiltonian (pH) input-state-output systems are defined as

Here, $J=-J^{\top}\in\mathbb{R}^{n\times n}$ is a structure matrix corresponding to energy-routing in the system, the matrix $0\leq R=R^{\top}\in\mathbb{R}^{n\times n}$ encodes the internal energy dissipation and $B\in\mathbb{R}^{n\times m}$ is an input matrix. The matrix $0\leq Q=Q^{\top}\in\mathbb{R}^{n\times n}$ captures the energy storage by means of the Hamiltonian $H(x):=\frac{1}{2}x^{\top}Qx$.

Due to the structure of the involved matrices, the trajectories of (pH) satisfy a power balance. More precisely, for all $t\geq 0$

If we integrate this equation, we obtain the corresponding energy balance, i.e., for all $t\geq 0$,

From these equations it can be seen that the supplied energy to the system is given by the passivity supply rate $u(t)^{\top}y(t)$. Regardless of the physical domain, e.g., mechanics, electrical engineering or fluid dynamics, cf. [26], the product of input and output will always have the physical unit of power, allowing for a coupling of port-Hamiltonian systems across different domains. If this coupling is performed in a power-preserving manner, the coupled system can be shown to be again port-Hamiltonian.

This section describes the modeling of adaptive structures, such as adaptive high-rise buildings. Structures comprising beams and links can be seen as flexible multi-body systems with infinitely many degrees of freedom. Such models include, e.g., Timoshenko or Euler-Bernoulli beam models and are given in the form of a partial differential equation. To obtain a model with finitely many degrees of freedom (DOFs), the finite element method is commonly applied. In this regard, a port-Hamiltonian model can be obtained by two means. First, the structural elements (links and beams here) are each modeled in port-Hamiltonian form, then spatially discretized using the partitioned finite element method (PFEM) from [2] and last assembled to a complex structure [29]. Second, a model in port-Hamiltonian form can be obtained by first deriving the linear second order equations of motion of the whole structure using the conventional finite element method (FEM) and then transforming it to port-Hamiltonian form.

Using the conventional FEM, the structures are discretized by finite beam or link elements, which are connected at nodes. Each node has a certain number of translational and rotational DOFs, depending on the configuration and used types of finite elements (e.g., links need only translational DOFs), the dimensionality (e.g.,two or three-dimensional) and constraints (e.g., at the support and coupling interfaces). The DOFs are collected in the vector $q(t)\in\mathbb{R}^{n_{\mathrm{dof}}}$ and the linear second order equations of motion are, for given initial DOFs $q_{0}\in\mathbb{R}^{n_{\mathrm{dof}}}$ and velocities $v_{0}\in\mathbb{R}^{n_{\mathrm{dof}}}$, given by

where $\dot{q}(t)$ are the nodal velocities, $\ddot{q}(t)$ are nodal accelerations and $M,K\in\mathbb{R}^{n_{\mathrm{dof}}\times n_{\mathrm{dof}}}$ are the positive definite mass and stiffness matrices, respectively.

Any kind of civil engineering structure exhibits some kind of damping. Rayleigh-damping as introduced in [24] is a common damping model in the literature. This is due to its numerically beneficial assumption that damping is a linear combination of mass and stiffness, $D=\alpha_{1}M+\alpha_{2}K$. For this choice, the undamped and the damped system exhibit identical eigenmodes. For the simulations in this paper we use $\alpha_{1}=0.05,\alpha_{2}=0.005$. For the proposed control method, a different damping model satisfying $D=D^{\top}\succ 0$ would also be feasible. The forces acting on the structure are collected in $f(t)\in\mathbb{R}^{n_{\mathrm{dof}}}$ and comprise external loads, such as wind loads and the forces applied by the actuators. In this work, no external loads are considered and a structure is excited by a non-zero initial condition, corresponding, e.g., to a previous displacement by wind.

The actuators are integrated in parallel, meaning that the actuator forces act in parallel to the hosting elements in opposite directions at the adjacent nodes. Thus $f(t)$ is rewritten as $f(t)=F_{\mathrm{u}}u(t)$, where $u(t)\in\mathbb{R}^{m}$ contains the actuator forces (the input) and the input matrix $F_{\mathrm{u}}\in\mathbb{R}^{{n_{\mathrm{dof}}}\times m}$ maps the input to nodal forces. For simplicity of notation, the time dependency and the initial conditions are omitted in the following.

The model equations (2.2) are transformed to port-Hamiltonian form (pH) using the DOFs and momenta as state variables $x=\left(\begin{smallmatrix}M\dot{q}\\ q\end{smallmatrix}\right)\in\mathbb{R}^{2n_{\mathrm{dof}}}$:

Thus, the energy

consists of kinetic and potential energy. Further, the collocated output according to (pH) is defined by

and has the physical unit of a velocity.

First, we inspect solutions of a linear quadratic problem given by

where $x^{0}\in\mathbb{R}^{n}$ is an initial state, $\mu>0$ is a regularization parameter and $\mathbb{U}\subset\mathbb{R}^{m}$ is a compact and convex control constraint set. In particular, $\mathbb{U}$ is bounded such that $L^{p}(0,T;\mathbb{U})\subset L^{\infty}(0,T;\mathbb{U})$ for all $1\leq p\leq\infty$.

In the remainder of this section, we state two theoretical results w.r.t. the long-time behavior of (3.1) in the form of a turnpike property. Such turnpike properties are crucial in the analysis of model predictive control schemes [7, 8]. We first consider the case of inactive control constraints, which, in view of Riccati theory is present if the penalization parameter $\mu>0$ is large w.r.t. the norm of the initial value.

Complementary to the case of inactive control constraints considered in Proposition 3.1, we now provide a result which also applies to the case of active control constraints.

We briefly provide an interpretation of the integral turnpike property of the previous proposition, i.e., Inequality (3.3). As the upper bound is independent of $T$ and as the integral is nonnegative, this means that if we increase the time horizon and thus the integration domain, the integrand on the left hand side has to be small for the majority of the time. This means, that such a bound provides us with a notion of smallness of the optimal state and control. Moreover, passing to the limit via $T\to\infty$ in (3.3) we obtain $x\in L^{2}(0,\infty;\mathbb{R}^{n})$ which implies the convergence $x(t)\stackrel{{\scriptstyle t\to\infty}}{{\to}}0$, cf. [33, Lemma 14.1].

Having provided stability results for the quadratic control penalization, we now introduce an alternative cost functional, which takes the port-Hamiltonian structure into account.

In order to include the pH-structure in the cost function and to minimize the supplied energy, cf. (2.1) we replace the quadratic control term in (3.1) by the supplied energy. Thus, for $x^{0}\in\mathbb{R}^{n}$, we consider the problem

For the case of energy-optimal state transitions, that is, in case $W=0$ and an additional terminal condition $x(T)=x_{T}\in\mathbb{R}^{n}$, this problem was analyzed in [6, 20, 5].

The integrand of the cost function of problem (3.5) is a particular case of the cost function of the LQR problem with cross-weighted costs as described in [13] with a certain choice of weighting matrices. In view of the the notation in [13], the cost function for such cross-weighted stage costs is given by $x(t)^{\top}Wx(t)+2x(t)^{\top}W^{1/2}Du(t)+u(t)^{\top}(P+D^{\top}D)u(t)$. The cost function integrand in (3.5) results using $D=\frac{1}{2}W^{-1/2}QB$ and $P=-D^{\top}D$. Here it is important to note that we do not have a quadratic control penalization, which leads, in general, to a singular optimal control problem. Invoking the port-Hamiltonian structure allows to resolve this issue and fully characterize the optimal state feedback, see [5].

Using the balance equation (2.1), the cost of (3.5) can be cast into the form

As $H(x(0))=H(x^{0})$, we will omit this term in the following, as it poses a mere offset of the cost functional. By means of the above reformulation, we see that the cost functional in (3.5) implicitly carries a terminal cost by means of the Hamiltonian. Moreover, a particular focus on the velocity components is intrinsically included by means of the dissipation matrix $RQ$ of (2.3).

Note that in (3.5) we did not include a penalization parameter $\mu>0$ for the control cost as for the quadratic cost in (3.1). The reason for this is simple: as can be observed in the numerical results in Section 4, we do not need to include a tuning of this cost functional to observe active damping. This is contrary to the quadratic control penalization in (3.1) where for too large penalization parameters no vibration damping beyond the natural dissipation is obtained.

First, similar to Propositions 3.1 and 3.2, we provide a result showing that energy-optimal solutions, i.e., solutions of (3.5) also stabilize the system. To this end, we denote the distance to a closed subset $M\subset\mathbb{R}^{n}$ by $\operatorname{dist}(x,M)=\min_{z\in M}\|x-z\|$.

For the quadratic control cost in (3.1), by standard Riccati theory and in the absence of control constraints, the control is given as a state feedback. For finite time horizons, this feedback can be obtained solving the differential Riccati equation, and on the infinite time horizon, by the algebraic Riccati equation, respectively. Due to the linear dependence of the cost on the control, such results can not be applied to (3.5). In the following, however, we show under a condition on the input matrix $B$, the dissipation matrix $RQ$ and the state weighting matrix $W$ in the cost that the control is given by means of state and adjoint state. This paves the way for future work considering state-feedback controls, as provided in [5] for state transition problems.

To derive optimality conditions, we define the optimality Hamiltonian

By Pontryagin’s maximum principle, cf., e.g., [12], optimal solutions $(x^{\star},u^{\star})\in W^{1,1}(0,T;\mathbb{R}^{n})\times L^{1}(0,T;\mathbb{U})$ of (3.5) satisfy

In view of [14, Rem. 6.9, p. 168], as we do not impose any terminal constraints, the multiplier $\lambda_{0}$ does not vanish such that we may normalize $\lambda_{0}=1$ in the following. Inserting the output equation $y=B^{\top}Qx$

and thus we obtain the optimality system

For components of the switching function

which do not vanish, the corresponding optimal control components are determined by means of a linear minimization problem and lie on the boundary of the compact control constraint set $\mathbb{U}$. However, if the $i$-th component of the switching function vanishes, the inclusion in the last line of (3.8) does not provide any information on the control component $u_{i}$, $i\in\{1,\ldots,m\}$. Thus, time intervals of positive measure where $s_{i}(t)\equiv 0$, $i\in\{1,\ldots,m\}$ are called singular arcs. On these singular arcs and on first sight, the optimality condition does not characterize the optimal control. The following result shows that under a certain condition linking the input matrix, the dissipation matrix, and the state weighting in the cost, this singularity can be resolved. More precisely, the optimal control on singular arcs can be expressed in terms of the state and the adjoint.

We conclude this section by a brief discussion of the assumption $\operatorname{im}B\cap\ker RQ\cap\ker W=\{0\}$ of Theorem 3.4. This condition is always satisfied (i.e., in particular for $W=0$) for fully damped (i.e., $D\succ 0$) second order systems of the form (2.3), as here already $\operatorname{im}B\cap\ker RQ=\{0\}$. This is a consequence of direct actuation only in the damped coordinates, i.e., the momenta. For $W=0$, it has been shown in [5] that the condition $\operatorname{im}B\cap\ker RQ=\{0\}$ is only sufficient and not necessary for characterizing the control by the state and the adjoint. More precisely, assuming full column rank of $B$, $\operatorname{im}B\cap\ker RQ=\{0\}$ holds if and only if the optimality system along the singular arc is a regular differential algebraic equation (DAE) with index three [5, Proposition 3.10]. We note that formula for the optimal control can be explicitly stated if the resulting DAE is regular, without a condition on the index [5].

First, we consider penalization of the full Hamiltonian, i.e., we solve OCPs (3.1) and (3.5) with $W=Q$, i.e.,

with penalization parameter $\mu>0$ and

The results in view of vibration damping are shown in Figure 3. We observe in the upper left plot that the port-Hamiltonian cost without tuning parameter outperforms all considered penalizations of the quadratic cost in view of the withdrawn energy. The corresponding numerical values are also displayed in Table 1. Further, in the second and third plot of Figure 3, it can be seen that the quadratic cost for $\mu=10^{-8}$ as well as the port-Hamiltonian cost introduce significant additional damping into the system. In contrast, the quadratic cost for $\mu=10^{-6}$ does not achieve a significant improvement compared to the uncontrolled cost. This reflects again an issue in the quadratic formulation: in order to obtain a good performance in view of vibration damping, a significant amount of tuning is necessary. A reason for this is that the squared control input term lacks physical interpretation. Contrary to this, in the port-Hamiltonian cost functional, both terms are intrinsic to the port-Hamiltonian structure and linked by the energy balance.

In the lower right plot of Figure 3 we also observe that for the quadratic cost and the choice $\mu=10^{-6}$ the control constraints are inactive over the full horizon. The reduced penalization $\mu=10^{-8}$ then leads to active control constraints on the initial part of the horizon, similar to the port-Hamiltonian cost functional. In particular this means that in the case of $\mu=10^{-6}$, a control obtained as a Riccati feedback is feasible, whereas this is no longer the case for $10^{-8}$.

Here, we consider an example with a weighting matrix for the state which only penalizes the Hamiltonian of the upmost level in the OCPs (3.1) and (3.5), i.e., its kinetic and potential energy. The results are depicted in Figure 4. First we see observe again, that the formulation with the port-Hamiltonian cost leads to a significantly higher withdrawn energy. In particular, we observe that for the smallest penalization parameter, the performance of the quadratic cost decreases again. This is contrary to the upper plot in Figure 3, where we observed a monotone behavior of the withdrawn energy in the penalization parameter. The precise values are depicted in Table 2.

Further, in the middle plot of Figure 4, we see that the quadratic cost with the depicted penalization parameters and the port-Hamiltonian cost induce sufficient damping in the system. We note that this is not true for higher damping parameters, where the higher control cost interferes with the goal of damping. We stress again, that such a tuning is not necessary for the port-Hamiltonian cost.