Introduction to "RORPack": A Python Software Library for Robust Output Regulation

TL;DR

RORPack is an open-source Python library that facilitates robust output regulation for linear PDE systems, including controller design, simulation, and visualization, with practical examples on heat and wave equations.

Contribution

This paper introduces RORPack, a comprehensive software library that automates controller design and simulation for robust output regulation of PDE systems.

Findings

Provides automated controller construction for PDEs

Includes simulation and visualization tools

Offers example cases on heat and wave equations

Abstract

This document contains the mathematical introduction to RORPack - a Python software library for robust output tracking and disturbance rejection for linear PDE systems. The RORPack library is open-source and freely available at https://github.com/lassipau/rorpack/ The package contains functionality for automated construction of robust internal model based controllers, simulation of the controlled systems, visualisation of the results, as well as a collection of example cases on robust output regulation of controlled heat and wave equations.

Equations

$$\dot{x}(t)$$

$$y(t)$$

$$\underset{t\to \infty }{\lim }\Vert y(t)-y_{\text{mbox}ref}(t){\Vert }_Y=0$$

$$y_{\text{mbox}ref}(t)$$

$$w_{\text{mbox}dist}(t)$$

$$\dot{z}(t)$$

$$u(t)$$

$$\frac{\partial x}{\partial t}(\xi ,t)$$

$$-\frac{\partial x}{\partial \xi }(0,t)$$

$$y(t)$$

$$\frac{\partial x}{\partial t}(\xi ,t)$$

$$-\frac{\partial x}{\partial \xi }(0,t)$$

$$y(t)$$

$$\frac{\partial x}{\partial t}(\xi ,t)$$

$$-\frac{\partial x}{\partial \xi }(0,t)$$

$$y(t)$$

$$b_1(\xi )$$

$$c_1(\xi )$$

$$x_t(\xi ,t)$$

$$\frac{\partial x}{\partial n}(\xi ,t){|}_{{\Gamma }_1}$$

$$y_1(t)$$

$$x_t(\xi ,t)$$

$$\frac{\partial x}{\partial n}(\xi ,t){|}_{{\Gamma }_1}$$

$$y(t)$$

$$x_{tt}(\xi ,t)$$

$$-\frac{\partial x}{\partial \xi }(0,t)$$

$$y(t)$$

$$x_{tt}(\xi ,t)$$

$$x(0,t)$$

$$x(\xi ,0)$$

$$y(t)$$

Full text

1. 1 General Description

2. 1.1 Controller Design and Simulation Workflow

3. 2 Introduction to Robust Output Regulation

4. 3 Using the Software

5. 4 Implemented Controller Types

6. 4.1 Comments on the Controller Parameters

7. 5 PDE Models of the Example Cases

8. 5.1 The 1D Heat Equations

9. 5.2 The 2D Heat Equations on Rectangular Domains

10. 5.3 The 1D Wave Equations

11. 6 Contributors

12. A External Links

1. General Description

RORPack is a Python library for controller design and simulation of robust output tracking and disturbance rejection for linear partial differential equation models. The package contains a number of complete examples on robust controller design and simulation of PDE models of the following types:

• •

one-dimensional diffusion equations with either boundary or distributed control and observation.

• •

two-dimensional heat equation on a rectangular domain.

• •

one-dimensional wave equations with either distributed or boundary control and observation.

New examples will also be added in the future versions of the package.

The purpose of this document is to give a general introduction to the background theory of robust output regulation for linear PDE systems, to describe the mathematical models that are included as the example cases, and to document the usage of the software package on a general level.

The purpose of RORPack is to serve as a tool to illustrate the theory of robust output regulation for distributed parameter systems and it should not (yet) be considered as a serious controller design software. The developers of the software do not take any responsibility and are not liable for any damage caused through use of this software. In addition, at its present stage, the library is not optimized in terms of numerical accuracy. Any helpful comments on potential improvements of the numerical aspects will be greatly appreciated!

The RORPack software is distributed as open source under the GNU General Public License version 3 (see LICENSE.txt for detailed license and copyright information) and can freely downloaded at the address

https://github.com/lassipau/rorpack/

This document is published in arXiv.org. The website of the project is located at the address (hosted by GitHub Pages)

https://lassipau.github.io/rorpack/

All comments and suggestions for improvements are welcome! These can be sent directly to lassi.paunonen@tuni.fi.

1.1. Controller Design and Simulation Workflow

Basic workflow for robust controller design and simulation for a given system is that the user creates a main file with the following parts:

1. (1)

Calling of user-defined Python routines that return a numerical approximation of the linear PDE model. 2. (2)

Defining the reference and disturbance signals to be considered. 3. (3)

Calling of a RORPack routine for construction of a robust controller that is suitable for the type of PDE model. This typically involves choosing appropriate parameters for the controller construction. 4. (4)

Calling of RORPack routines for construction and simulation of the closed-loop system. The routines return numerical data describing the output, the error, the control signal, and the closed-loop state. 5. (5)

Calling of RORPack routines for visualising the behaviour of the output and the output error, as well as user-defined routines for illustrating the behaviour of the state of the controlled PDE system (multidimensional plots or animations).

The example cases included in RORPack are built around main files that follow the same structure, and the users are encouraged to implement their own simulations by modifying the example codes.

It is also possible to combine other Python numerical software packages in the study of robust controller design for PDE models in order to employ ready-made numerical approximations or additional numerical methods in the computation of required controller parameters. This approach is illustrated in one of the PDE example cases “heat_1d_2” described in Section 5.1 (Case 2). In this example the transfer function of the controlled PDE and other parameters used in the controller construction are computed with the Chebfun package [11, 41] (https://chebfun.org/) which provides flexible tools for solution of differential equations with extreme accuracy using spectral methods.

2. Introduction to Robust Output Regulation

In this section we give a general introduction to the mathematical theory of robust output regulation. The purpose of the RORPack package is to illustrate controller design for linear distributed parameter systems of the form

[TABLE]

on a Banach or a Hilbert space $X$. Controlled linear PDE models describing diffusion-convection phenomena, waves and vibrations and elastic deformations can be written in this form with a suitable differential operator $A$ [7, 42, 23]. In our main control problem the goal is to design a dynamic error feedback controller in such a way that the output $y(t)$ of the system converges to a given reference signal $y_{\mbox{\scriptsize{ref}}}(t)$ despite the external disturbance signal $w_{\mbox{\scriptsize{dist}}}(t)$. In addition, the controller needs to be robust in the sense that it achieves the output tracking and disturbance rejection even if the parameters $(A,B,B_{d},C,D,D_{d})$ are perturbed or contain small uncertainties.

Our main emphasis is on robust output regulation of diffusion-convection equations, wave equations, and beam and plate equations. However, the RORPack package can also be used for construction of controllers for finite-dimensional systems with given matrices $(A,B,B_{d},C,D,D_{d})$.

The full robust output regulation problem is defined in the following way.

The Robust Output Regulation Problem.* Given a reference signal $y_{\mbox{\scriptsize{ref}}}(t)$, design a dynamic error feedback controller such that the output $y(t)$ of the system converges to the reference signal asymptotically, i.e.,*

[TABLE]

*despite the disturbance signal $w_{\mbox{\scriptsize{dist}}}(t)$. Moreover, the controller is required to be robust in the sense that it achieves the convergence of the output (2.2) even under small uncertainties and changes in the parameters $(A,B,B_{d},C,$ $D,D_{d})$ of the system. *

The reference signal $y_{\mbox{\scriptsize{ref}}}(t)$ and the disturbance signals $w_{\mbox{\scriptsize{dist}}}(t)$ are assumed to be of the form

[TABLE]

for some known frequencies $\{\omega_{k}\}_{k=0}^{q}\subset{\mathbb{R}}$ with $0=\omega_{0}<\omega_{1}<\ldots<\omega_{q}$ and unknown amplitudes $\{a_{k}^{j}\}_{k,j},\{b_{k}^{j}\}_{k,j}\subset{\mathbb{R}}$ (some of which may zero).

Dynamic feedback is essential for achieving robust output regulation, and the control problem can indeed be solved with a dynamic error feedback controller (see Figure 1). The classical internal model principle gives a characterization for the controllers that solve the robust output regulation problem. This fundamental result was first introduced for finite-dimensional linear systems in the 1970’s by Francis and Wonham [12] and Davison [8] (see [17, Ch. 1] for an excellent overview). The internal model principle was later extended for infinite-dimensional linear systems by the Systems Theory Research Group at Tampere, Finland in the references [20, 34, 29, 36, 30]111Our focus is on linear systems, but there is also an extensive literature on the internal model principle for nonlinear systems, see, e.g., [17, 4, 5] and references therein..

The internal model principle [34, Thm. 6.9] is also the most important tool in designing controllers for robust output regulation. The result states that a controller solves the robust output regulation problem if and only if the following conditions are satisfied.

• •

The error feedback controller incorporates “an internal model” of the signals $y_{\mbox{\scriptsize{ref}}}(t)$ and $w_{\mbox{\scriptsize{dist}}}(t)$ in (2.3).

• •

The closed-loop system is exponentially or strongly stable.

In controller design, the internal model property can be guaranteed by choosing a suitable structure for the controller. The rest of parameters are subsequently chosen so that the closed-loop system becomes stable.

The detailed description of the theory of robust output regulation problem and the internal model principle can be found in the references listed below. The main emphasis in the list is (shamelessly) on the publications by the Systems Theory Research Group at Tampere University, Finland.

• •

[38, 43, 13, 24, 14, 39, 25, 22, 21, 10, 2, 15, 16, 35, 44] (including PI-control for PDE models): Robust controller design in various forms for infinite-dimensional linear systems.

• •

[34] and [36]: The Internal Model Principle for infinite-dimensional linear systems with bounded and unbounded, respectively, input and output operators $B$ and $C$.

• •

[30, 31]: Robust controller design for regular linear systems.

• •

[19, 18]: Robust controller design for port-Hamiltonian systems and other boundary controlled partial differential equations.

• •

[39, 32, 33]: Robust controller design impedance passive systems.

• •

[28] for robust finite-dimensional low-order controller design for parabolic systems using Galerkin approximations and model reduction.

• •

References on the output regulation without the robustness requirement: [37, 40, 6, 9, 27, 45]

The considered controllers are of the form

[TABLE]

Here $y(t)-y_{\mbox{\scriptsize{ref}}}(t)$ is the regulation error. The construction of the robust controllers are based on the references [30, 32] with certain modifications and improvements. In particular, the internal models are defined in their “real forms”, making the controller real whenever the parameters of the plant are real. The same controllers also achieve robust output tracking also for reference and disturbance signals with complex coefficients $\{a_{k}^{j}\}_{k,j},\{b_{k}^{j}\}_{k,j}\subset{\mathbb{C}}$.

The constructions of the robust controllers use the knowledge of the frequencies $\{\omega_{k}\}_{k}$ of the reference and disturbance signals, the number of outputs $p:=\dim Y$ of the system (2.1), and certain knowledge of the system. In particular, the “minimal controllers” require knowledge of the values $P(i\omega_{k})$ of the transfer function $P(\lambda)=C_{\Lambda}R(\lambda,A)B+D$ of the system at the complex frequencies $\{i\omega_{k}\}_{k}\subset i{\mathbb{R}}$ of the reference and disturbance signals (2.3). The other controller structures also use knowledge of the parameters $(A,B,C,D)$ of the system as they involve designing an observer for (2.1), and require the user to provide stabilizing state feedback and output injection operators. See Section 4 for more information on the controller-specific requirements.

3. Using the Software

The rorpack Python library can be installed for Python 2 or Python 3 as instructed in the README.md file included in the software package (this mainly involves downloading the source codes and typing ’pip install .’ or ’pip3 install .’ in the main directory). The subdirectory examples/ contains the included PDE examples and simulation files (documented in detail in Section 5).

The following commented example file explains the typical structure and workflow of the controller construction and simulation with RORPack. The considered example cases included in the file heat_1d_3.py and documented in Section 5.1. Due to the properties of the Python language, the constructor routines used in the main simulation are written at the beginning of the file.

Contents of the file heat_1d_3.py: The file begins with comment lines and loading of the necessary parts of the rorpack library as well as other Python packages.

⬇

”’

Heat equation on the interval [0,1] with Neumann boundary control

and Dirichlet boundary observation. Approximation with a Finite

Difference scheme.

Neumann boundary disturbance at x=0, two distributed controls and

two distributed measurements regulated outputs. The controls act

on the intervals ’IB1’ and ’IB2’ (Default ’IB1’ = [0.3,0.4] and

’IB2’ = [0.6,0.7]) and the measurements are the average

temperatures on the intervals ’IC1’ and ’IC2’ (Default

’IC1’ = [0.1,0.2] and ’IC2’ = [0.8,0.9]).

’”

import numpy as np

from rorpack.system import LinearSystem

from rorpack.controller import *

from rorpack.closed_loop_system import ClosedLoopSystem

from rorpack.plotting import *

from laplacian import diffusion_op_1d

The next part introduces a constructor routine to define the Finite Difference approximation of the heat equation and the input and output operators. The parameter $N$ is the size of the approximation and cfun is a function describing the spatially varying thermal diffusivity of the material.

⬇

def construct_heat_1d_3(N, cfun, IB1, IB2, IC1, IC2):

spgrid = np.linspace(0, 1, N+1)

plt.plot(spgrid,cfun(spgrid))

plt.title(’The thermal diffusivity $c(x)$ of the material’)

plt.tight_layout()

plt.grid(True)

plt.show()

h = spgrid[1]-spgrid[0]

DiffOp, spgrid = diffusion_op_1d(spgrid, cfun, ’ND’)

A = DiffOp

B1 = 1/(IB1[1] - IB1[0])*np.logical_and(spgrid >= IB1[0],

        spgrid <= IB1[1])

B2 = 1/(IB2[1] - IB2[0])*np.logical_and(spgrid >= IB2[0],

        spgrid <= IB2[1])

B = np.stack((B1, B2), axis=1)

C1 = h/(IC1[1] - IC1[0])*np.logical_and(spgrid >= IC1[0],

        spgrid <= IC1[1])

C2 = h/(IC2[1] - IC2[0])*np.logical_and(spgrid >= IC2[0],

        spgrid <= IC2[1])

C = np.stack((C1, C2))

D = np.zeros((2, 2))

Bd = np.bmat([[np.atleast_2d(-2/h)], [np.zeros((N-1, 1))]])

return LinearSystem(A, B, C, D, Bd, np.zeros((2,1))), spgrid

The next part defines the parameters of the system and constructs the system $(A,B,B_{d},C,D)$ as an object of the class LinearSystem of RORPack.

⬇

Parameters for this example.

N = 100

The spatially varying thermal diffusivity of the material

cfun = lambda x: np.ones(np.atleast_1d(x).shape)

cfun = lambda x: 1+x

cfun = lambda x: 1-2x(1-2*x)

cfun = lambda x: 1+.5np.cos(5/2np.pi*x)

Note: Lower diffusivity is difficult for the Low-Gain

and Passive controllers

cfun = lambda x: .2-.4x(1-x)

Regions of inputs and outputs

IB1 = np.array([0.3, 0.4])

IB2 = np.array([0.6, 0.7])

IC1 = np.array([0.1, 0.2])

IC2 = np.array([0.8, 0.9])

Length of the simulation

t_begin = 0

t_end = 8

t_points = 300

Construct the system.

sys, spgrid = construct_heat_1d_3(N, cfun, IB1, IB2, IC1, IC2)

The next part defines the reference signal $y_{\mbox{\scriptsize{ref}}}(t)$ (in yref) and the disturbance signal $w_{\mbox{\scriptsize{dist}}}(t)$ (in wdist) and lists the (real) frequencies $\{\omega_{k}\}_{k=1}^{q}$ in the variable freqsReal. Alternative reference and disturbance signals are commented out in the code for further simulation experiments.

⬇

Define the reference and disturbance signals, and list the

required frequencies in ’freqsReal’

Case 1:

yref = lambda t: np.stack((np.sin(2t), 2np.cos(3*t)))

wdist = lambda t: np.sin(6*t)

Case 2:

yref = lambda t: np.ones((2,np.atleast_1d(t).shape[0]))

wdist = lambda t: np.zeros(np.atleast_1d(t).shape)

freqsReal = np.array([0, 1, 2, 3, 6])

The next part constructs the chosen controller structure and the closed-loop system as objects of the RORPack classes. Alternative controller structures are commented out in the code for easy comparison of controller performances.

⬇

Construct the controller and the closed loop system.

Controller choices, Low-gain robust controller

Requires the transfer function values P(i*w_k)

epsgainrange = np.array([0.3,0.6])

Pvals = np.array(list(map(sys.P, 1j * freqsReal)))

contr = LowGainRC(sys, freqsReal, epsgainrange, Pvals)

Dual observer-based controller

Requires stabilizing operators K and L1

and the transfer function values P_L(i*w_k)

K = -sys.B.conj().T

L1 = -10*sys.C.conj().T

PLvals = np.array(list(map(lambda freq: sys.P_L(freq, L1),

1j * freqsReal)))

IMstabmargin = 0.5

IMstabmethod = ’LQR’

contr = DualObserverBasedRC(sys, freqsReal, PLvals, K, L1,

IMstabmargin, IMstabmethod)

Observer-based controller

Requires stabilizing operators K21 and L

and the transfer function values P_K(i*w_k)

K21 = -sys.B.conj().T

L = -10*sys.C.conj().T

PKvals = np.array(list(map(lambda freq: sys.P_K(freq, K21),

1j * freqsReal)))

IMstabmargin = 0.5

IMstabmethod = ’LQR’

contr = ObserverBasedRC(sys, freqsReal, PKvals, K21, L,

    IMstabmargin, IMstabmethod)

Construct the closed-loop system

clsys = ClosedLoopSystem(sys, contr)

The next part simulates the behaviour of the closed-loop system. The initial state of the system is chosen in the variable x0, and the initial state of the controller is chosen by default to be zero.

⬇

Simulate the system, define the initial state x_0

x0fun = lambda x: np.zeros(np.size(x))

x0fun = lambda x: 0.5 * (1 + np.cos(np.pi * (1 - x)))

x0fun = lambda x: 3*(1-x)+x

x0fun = lambda x: 1/2x**2(3-2*x)-1

x0fun = lambda x: 1/2x**2(3-2*x)-0.5

x0fun = lambda x: 1*(1-x)2(3-2(1-x))-1

x0fun = lambda x: 0.5*(1-x)2(3-2(1-x))-0.5

x0fun = lambda x: 0.25*(x3-1.5*x2)-0.25

x0fun = lambda x: 0.2x**2(3-2*x)-0.5

x0 = x0fun(spgrid)

z0 is chosen to be zero by default

z0 = np.zeros(contr.G1.shape[0])

xe0 = np.concatenate((x0, z0))

tgrid = np.linspace(t_begin, t_end, t_points)

sol,output,error,control,t = clsys.simulate(xe0,tgrid, yref, wdist)

print(’Simulation took %f seconds’ % t)

Finally, the results of the simulation are plotted in separate figures and the behaviour of the controlled state is animated using the user-defined routines.

⬇

Plot the output and the error, and animate the behaviour

of the controlled state.

plot_output(tgrid, output, yref, ’subplot’, ’default’)

plot_error_norm(tgrid, error)

plot_control(tgrid, control)

In plotting and animating the state, fill in the homogeneous

Dirichlet boundary condition at x=1

sys_state = np.vstack((sol.y[0:N],np.zeros((1,np.size(tgrid)))))

spgrid = np.concatenate((spgrid,np.atleast_1d(1)))

plot_1d_surface(tgrid, spgrid, sys_state, colormap=cm.plasma)

animate_1d_results(spgrid, sys_state, tgrid)

4. Implemented Controller Types

In this section we list the concrete controllers implemented in RORPack. The documentation of the code includes additional information on the usage of the construction routines. The controllers are the following:

• •

The “minimal robust controller” (including only the internal model), based on references [14], [30, Sec. IV]. Stabilization of the closed-loop system is based on selection of a suitable low-gain parameter $\varepsilon>0$.

Calling sequence for the construction:

LowGainRC(sys,freqsReal,epsgain,Pvals)

where sys contains the parameters of the plant, freqs contains the frequencies $\{\omega_{k}\}_{k=0}^{q}$ of the signals $y_{\mbox{\scriptsize{ref}}}(t)$ and $w_{\mbox{\scriptsize{dist}}}(t)$ in (2.3), epsgain is the value of the low-gain parameter $\varepsilon>0$. Finally, Pvals is a $(q+1)\times p\times m$ array containing the values $P(i\omega_{k})=C_{\Lambda}R(i\omega_{k},A)B+D\in{\mathbb{C}}^{p\times m}$ of the transfer function of (2.1) at the complex frequencies $\{i\omega_{k}\}_{k=0}^{q}$ of the reference and disturbance signals in (2.3). The parameter epsgain can alternatively be a vector of length 2 providing minimal and maximal values for $\varepsilon$. The controller construction has a naive functionality for finding an $\varepsilon$ to optimize stability margin of the numerically approximated closed-loop system (simply by starting from the minimal value and increasing $\varepsilon$ in steps).

• •

The “observer-based robust controller”, based on references [15, Sec. 7], [30, Sec. VI], [31, Sec. 5]. The closed-loop stability is achieved using an observer for the state of the system (2.1).

Calling sequence for the construction:

ObserverBasedRC(sys,freqsReal,PKvals,K21,L,

a                        IMstabmargin,IMstabmethod,CKRKvals)

where sys contains the parameters of the plant, freqsReal contains the frequencies $\{\omega_{k}\}_{k=0}^{q}$ of the signals $y_{\mbox{\scriptsize{ref}}}(t)$ and $w_{\mbox{\scriptsize{dist}}}(t)$ in (2.3). The parameters K21 and L describe operators $K_{21}$ and $L$, respectively, such that $A+BK_{21}^{\Lambda}$ and $A+LC_{\Lambda}$ generate exponentially stable semigroups. The parameter PKvals is a $(q+1)\times p\times m$ array containing the values $P_{K}(i\omega_{k})=(C_{\Lambda}+DK_{21}^{\Lambda})R(i\omega_{k},A+BK_{21}^{\Lambda})B+D\in{\mathbb{C}}^{p\times m}$ of the stabilized transfer function of (2.1) at the complex frequencies $\{i\omega_{k}\}_{k=0}^{q}$ of the reference and disturbance signals in (2.3).

Instead of the approach used in [30], the “internal model”, i.e., the pair $(G_{1},B_{1})$, in the controller is stabilized using either LQR-based design (IMstabmethod = ’LQR’, by default) or pole placement (IMstabmethod = ’poleplacement’) with a stability margin IMstabmargin (default = 0.5). Note that the variable IMstabmargin only determines the stability margin of the internal model, and the stability margin of the closed-loop system also depends on the stability margins of the semigroups generated by $A+BK_{2}^{\Lambda}$ and $A+LC_{\Lambda}$. Finally, CKRKvals (optional) is a $(q+1)\times p\times N$ array containing elements $(C_{\Lambda}+DK_{21}^{\Lambda})R(i\omega_{k},A+BK_{21}^{\Lambda})\in{\mathbb{C}}^{p\times N}$ for $k=\{0,\ldots,q\}$. If this parameter is not given, the routine computes these values based on the same numerical approximation as the one used in the simulation.

• •

The “dual observer-based robust controller”, based on references [30, Sec. V], [31, Sec. 4]. The closed-loop stability is achieved using a complementary controller structure that coincides with observer-based stabilization of the dual of the closed-loop system.

Calling sequence for the construction:

DualObserverBasedRC(sys,freqsReal,PLvals,K2,L1,

a                        IMstabmargin,IMstabmethod,RLBLvals)

where sys contains the parameters of the plant, freqs contains the frequencies $\{\omega_{k}\}_{k=0}^{q}$ of the signals $y_{\mbox{\scriptsize{ref}}}(t)$ and $w_{\mbox{\scriptsize{dist}}}(t)$ in (2.3). The parameters K2 and L1 describe operators $K_{2}$ and $L_{1}$, respectively, such that $A+BK_{2}^{\Lambda}$ and $A+L_{1}C_{\Lambda}$ generate exponentially stable semigroups. The parameter PLvals is a $(q+1)\times p\times m$ array containing the values $P_{L}(i\omega_{k})=C_{\Lambda}R(i\omega_{k},A+L_{1}C_{\Lambda})(B+L_{1}D)+D\in{\mathbb{C}}^{p\times m}$ of the stabilized transfer function of (2.1) at the complex frequencies $\{i\omega_{k}\}_{k=0}^{q}$ of the reference and disturbance signals in (2.3).

Instead of the approach in [30], the “internal model”, i.e., the pair $(C_{1},G_{1})$, in the controller is stabilized using either LQR-based design (IMstabmethod = ’LQR’) or pole placement (IMstabmethod = ’poleplacement’) with a predefined stability margin IMstabmargin. Note that the variable IMstabmargin only determines the stability margin of the internal model, and the stability margin of the closed-loop system also depends on the stability margins of the semigroups generated by $A+BK_{2}^{\Lambda}$ and $A+L_{1}C_{\Lambda}$.

Finally, RLBLvals (optional) is a $(q+1)\times N\times m$ array containing elements $R(i\omega_{k},A+L_{1}C_{\Lambda})(B+L_{1}D)\in{\mathbb{C}}^{N\times m}$ for $k=\{0,\ldots,q\}$. If this parameter is not given, the routine computes these values based on the same numerical approximation as the one used in the simulation.

• •

The “passive minimal controller” based on [39, Thm. 1.2], [32, Sec. 5.1], [33]. For a strongly stable impedance passive system the closed-loop stability can be achieved using passive controller design and a power preserving interconnection between the passive control system and the controller.

Calling sequence for the construction:

PassiveRC(freqsReal,dim_Y,epsgain,sys)

where freqsReal contains the (real) frequencies $\{\omega_{k}\}_{k=0}^{q}$ of the signals $y_{\mbox{\scriptsize{ref}}}(t)$ and $w_{\mbox{\scriptsize{dist}}}(t)$ in (2.3), dim_Y is the dimension of the output space, epsgain is a parameter $\varepsilon>0$ that controls the norm of $K$, and sys contains the parameters of the plant. In addition, Pvals is a $(q+1)\times p\times m$ array containing the values $P(i\omega_{k})=C_{\Lambda}R(i\omega_{k},A)B+D\in{\mathbb{C}}^{p\times m}$ of the transfer function of (2.1) at the complex frequencies $\{i\omega_{k}\}_{k=0}^{q}$ of the reference and disturbance signals in (2.3). The parameter Pvals is optional and it is currently not used in the controller design. In future versions of the implementation this information may be used to improve the performance of the controller.

The parameter epsgain can alternatively be a vector of length 2 providing minimal and maximal values for $\varepsilon$. The controller construction has a naive functionality for finding an $\varepsilon$ to optimize stability margin of the numerically approximated closed-loop system (simply by starting from the minimal value and increasing $\varepsilon$ in steps).

It should be noted that the controller construction routine does not test passivity, and for a non-passive system the resulting closed-loop system will typically be unstable.

4.1. Comments on the Controller Parameters

In this section we make some remarks on the choices of the controller parameters.

**The gain parameter $\varepsilon>0$ in **LowGainRC. The theory in [14, 39] guarantees that for a stable system $(A,B,C,D)$ there exists $\varepsilon^{\ast}>0$ such that for any $0<\varepsilon<\varepsilon^{\ast}$ the closed-loop system is exponentially stable. The stability margin of the closed-loop system (which directly determines the convergence rate of the regulation error $e(t)$) can be optimized with a suitable choice of $\varepsilon>0$. Finding such an optimal value of $\varepsilon>0$ for a PDE system is a challenging task, but in numerical simulations one can use a naive approach of tracking the spectrum of the finite-dimensional closed-loop system matrix $A_{e}(\varepsilon)$. This is the approach taken in the example cases in RORPack, though it should be noted that for general PDE systems there is no guarantee that the value $\varepsilon>0$ obtained this way would optimize the closed-loop system for the original PDE system. In the example cases we are mainly interested in finding an $\varepsilon>0$ which achieves a reasonable rate of convergence rate of the error $\lVert e(t)\rVert$.

The gain in the passive robust controller PassiveRC can analogously be adjusted with a choice of a parameter $\varepsilon>0$ to optimize the closed-loop stability margin.

The transfer function values. The controllers make use of the values $P(i\omega_{k})$, $P_{K}(i\omega_{k})$ and $P_{L}(i\omega_{k})$ of the transfer functions associated to the original or stabilized versions of the PDE system. The RORPack class LinearSystem has the necessary routines for computing these values based on the matrices A, B, C, and D. However, using the same approximation of the PDE for both controller design and simulation corresponds to essentially controlling the approximation as a finite-dimensional system. If possible, to avoid unrealistically positive results, it is therefore better to use two different approximations for controller construction and simulation. This is especially the case if the validity of the approximation for the computation of the parameters can not be guaranteed with absolute certainty. However, there are of course cases, such Galerkin approximations of parabolic systems [26], where the values of the controller parameters can be shown to converge with the order of the approximation, and above concerns are unnecessary. The same comments apply to the parameters CKRKvals and RLBLvals.

For some special PDEs, such as 1D heat or wave equations with constant coefficients, the values $P(i\omega_{k})$ may have explicit expressions, but these are very limited special cases. For PDEs with spatially varying parameters there are powerful computational methods that can be used to determine $P(i\omega_{k})$, such as the Chebfun package (https://chebfun.org/) employed in heat_1d_2.py in Section 5.1 (“Case 2”).

5. PDE Models of the Example Cases

In this section we introduce the PDE models considered in the example cases and review their fundamental properties.

The considered reference and disturbance signals are in each case combinations of trigonometric functions of the form (2.3) with a given set of frequencies. The precise choices of the signals can be seen from the main files of the examples. Similarly, the chosen initial states are visible from the source files, and in several files alternative initial states are provided (these can be used by uncommenting the corresponding lines of code).

5.1. The 1D Heat Equations

This collection of examples consider 1D heat equations with spatially varying thermal diffusivity on $\Omega=[0,1]$ with different configurations of control inputs and measured outputs. The first two cases consider control, disturbance, and observation on the boundary, and the final example considers a system with two distributed inputs and outputs. In general, any combination of the above types of inputs and outputs is possible. The main property from the point of view of robust output regulation is whether or not the uncontrolled system is exponentially stable (minimal low-gain controller can be used) or impedance passive (the passive robust controller can be used).

In each of the cases, the semidiscretization of the PDE is completed using Finite Differences.

Case 1: Neumann boundary input at $\xi=0$, disturbance and output at $\xi=1$

Main file name: heat_1d_1.py

The model on $\Omega=[0,1]$ is

[TABLE]

where $c(\cdot)\geq c_{0}>0$ is the spatially varying thermal diffusivity of the material. The uncontrolled system is unstable due to the eigenvalue $\lambda=0$. The example uses the “Observer-Based Robust Controller” and “Dual Observer-Based Robust Controller” to achieve robust output tracking and disturbance rejection.

Case 2: Input, output, and disturbance at $\xi=0$, Dirichlet at $\xi=1$

Main file name: heat_1d_2.py

The model on $\Omega=[0,1]$ is

[TABLE]

where $c(\cdot)\geq c_{0}>0$ is the spatially varying thermal diffusivity of the material. The uncontrolled system is exponentially stable due to the homogeneous Dirichlet boundary condition at $\xi=1$. The system is also impedance passive since the control input and measured output are collocated, and because of this the robust output regulation problem can be solved using the Passive Robust Controller.

In this example we use the Chebfun** package** (https://chebfun.org/) [11, 41] in computing the values $P(i\omega_{k})$ of the transfer function as well as other parameters required in the controller construction. The Chebfun package utilizes spectral methods and provides powerful and easy-to-use tools for the solution of (especially 1D) boundary value problems with accuracies close to machine precision. Because of this, Chebfun has great potential in controller design for PDEs and it is especially perfectly suited for robust output regulation of this class of systems. The only drawback from RORPack’s perspective is that at the moment Chebfun is only implemented in Matlab222Though partial Python implementations exist, the most important functionality for solving BVPs can only be found in the original Matlab version.. However, using the Python-to-Matlab interface library “matlab.engine” included in Matlab (since R2014b), it is possible to call Matlab scripts and functions directly from Python. This approach is used by default in the current PDE example case (this requires installed versions of Matlab, Chebfun package, and the separate installation matlab.engine Python package, see https://se.mathworks.com/help/matlab/matlab-engine-for-python.html for details). The downside of this approach is that the startup time of the Matlab engine in Python can be extremely slow, and therefore the computations require a relatively long time. However, we get the benefit of very accurate values of the controlled PDE system to be used in the controller design.

Alternate computations of the controller parameters using the Finite Difference approximation (the one used in the main simulation) are commented out in the example code, and can be uncommented to run the heat_1d_2 example case without the Matlab interface.

Case 3: Distributed input and output, boundary disturbance at $\xi=0$, Dirichlet at $\xi=1$

Main file name: heat_1d_3.py

The model on $\Omega=[0,1]$ is

[TABLE]

where $c(\cdot)\geq c_{0}>0$ is the spatially varying thermal diffusivity of the material and

[TABLE]

Here $\chi_{[a,b]}(\cdot)$ is the characteristic function on the interval $[a,b]\subset[0,1]$, and thus the control inputs act on the intervals $[0.3,0.4]$ and $[0.6,0.7]$, and the outputs measure the average temperatures on the intervals $[0.1,0.2]$ and $[0.8,0.9]$. The uncontrolled system is stable because of the homogeneous Dirichlet boundary condition at $\xi=1$.

Figure 2 shows example results of the simulations including plots of the outputs and reference signals, norm of the regulation error, computed control signals, and the state of the controlled system as a function of $\xi$ and $t$.

5.2. The 2D Heat Equations on Rectangular Domains

These examples consider two-dimensional heat equations on $\Omega=[0,1]\times[0,1]$ with boundary input and output and boundary disturbances. The inputs and outputs act in an averaged sense on the boundaries.

The reference and disturbance signals are again of the form (2.3) and can be seen from the main files of the simulations.

Case 1: Collocated input and output

Main file name: heat_2d_1.py

This case is considered in [30, Sec. VII]. The system on $\Omega=[0,1]\times[0,1]$ with two inputs $u(t)=(u_{1}(t),u_{2}(t))^{T}$ and two outputs $y(t)=(y_{1}(t),y_{2}(t))^{T}$ is determined by

[TABLE]

Here the parts $\Gamma_{0}$, $\Gamma_{1}$, and $\Gamma_{2}$ of the boundary $\partial\Omega$ are defined so that $\Gamma_{1}=\{\,\xi=(\xi_{1},0)\mid 0\leq\xi_{1}\leq 1/2\,\}$, $\Gamma_{2}=\{\,\xi=(\xi_{1},1)\mid 1/2\leq\xi_{1}\leq 1\,\}$, $\Gamma_{0}=\partial\Omega\setminus(\Gamma_{1}\cup\Gamma_{2})$. By [3, Cor. 2] the heat equation defines a regular linear system with feedthrough $D=0$. The system is also impedance passive. The uncontrolled system is unstable due to the eigenvalue $\lambda=0$, but it can be stabilized exponentially with negative output feedback $u(t)=-\kappa y(t)$ for any $\kappa>0$. In this example we assume the system is pre-stabilized with $\kappa=1$.

In the simulations, the system is approximated using the eigenmodes of the Laplacian. The system can be controlled either with the low-gain minimal robust controller or the passive controller (with pre-stabilizing negative output feedback), or with one of the observer-based robust controllers.

Case 2: Non-collocated input and output

Main file name: heat_2d_2.py

A similar case (but with collocated inputs and outputs) was considered in [32, Sec. 6.3]. The system on $\Omega=[0,1]\times[0,1]$ is

[TABLE]

Here the parts $\Gamma_{0}$, $\Gamma_{1}$, $\Gamma_{2}$, and $\Gamma_{3}$ of the boundary $\partial\Omega$ are defined so that $\Gamma_{1}=\{\,\xi=(0,\xi_{2})\mid 0\leq\xi_{2}\leq 1\,\}$, $\Gamma_{2}=\{\,\xi=(\xi_{1},1)\mid 0\leq\xi_{1}\leq 1\,\}$, $\Gamma_{3}=\{\,\xi=(\xi_{1},0)\mid 0\leq\xi_{1}\leq 1/2\,\}$, $\Gamma_{0}=\partial\Omega\setminus(\Gamma_{1}\cup\Gamma_{2}\cup\Gamma_{3})$. By [3, Cor. 2] the heat equation defines a regular linear system with feedthrough $D=0$. The uncontrolled system is unstable due to the eigenvalue $\lambda=0$.

In the simulations, the system is approximated using Finite Differences with a uniform grid.

5.3. The 1D Wave Equations

In these examples we consider one-dimensional undamped and damped wave equations with control and observation at the boundaries and inside the domain.

Case 1: Non-collocated boundary input and output, disturbance near the output

Main file name: wave_1d_1.py

The model on $\Omega=[0,1]$ is

[TABLE]

The uncontrolled system is unstable due to lack of damping. It also cannot be stabilized with output feedback due to the non-collocated configuration of the inputs and outputs. The goal is to achieve tracking of the output $y(t)$. The additional measured output $y_{m}(t)$ is required to achieve closed-loop stability due to the fact that the system is not exponentially detectable with output $y(t)$ (because of the unobservability of the eigenvalue $\lambda=0$).

In the example case, the second output $y_{m}(t)$ is used to prestabilize eigenvalue $\lambda=0$ of the system with preliminary output feedback of the form $u(t)=-\kappa_{m}y_{m}(t)+\tilde{u}(t)$ (where $\tilde{u}(t)$ is the new input of the system). There are also other (and better) ways to handle the situation, and these will be implemented in later versions of this example. The pairs $(A,B)$ and $(C,A)$ are stabilized using collocated designs, i.e., we choose $K=-\kappa B^{\ast}$ and $L=-\ell C^{\ast}$ to stabilize the pairs $A+BK$ and $A+LC$. The main file contains the pre-stabilization gain $K_{m}>0$ and the gains $\kappa,\ell>0$ as design parameters.

The example considers output tracking of the velocity $x_{t}(0,t)$ to arbitrary 2-periodic reference signals. This is achieved by including frequencies of the form $k\pi$ for $k\in\{1,\ldots,q\}$ in the internal model. The reference signal is defined by defining its profile over one period in the variable ’yref1per’. Note that it is not necessary to compute the amplitudes of the frequency components of $y_{\mbox{\scriptsize{ref}}}(t)$ (which would be equivalent to finding the Fourier series expansion of the reference signal).

In the simulations the system is approximated using the orthonormal eigenfunctions of the undamped wave operator. Figure 3 shows example results of the simulation for two different periodic reference signals — a nonsmooth triangle signal and signal consisting of semicircles.

Case 2: Collocated distributed input and output

Main file name: wave_1d_2.py

The model on $\Omega=[0,1]$ is

[TABLE]

for some $b(\cdot),b_{d}(\cdot)\in L^{2}(0,1;{\mathbb{R}})$. The uncontrolled system is unstable due to the lack of damping. Since the input and output are distributed, it is also not exponentially stabilizable or detectable, but instead it is only strongly (or polynomially [1]) stabilizable provided that $\langle b(\cdot),\sin(k\pi\cdot)\rangle_{L^{2}}\neq 0$ for all $k\in{\mathbb{N}}$. Because of this, it is not possible to use the low-gain controller, and the two observer-based controller designs are not guaranteed to achieve closed-loop stability. However, since the system is impedance passive, the Passive Robust Controller can be used in achieving robust output regulation even in the absence of exponential stability as shown in [32, Sec. 5.1]. Due to the sub-exponential closed-loop stability, the regulation error does not converge with any uniform convergence rate, but instead the rate depends on the initial state of the system.

In the simulations the system is approximated using the orthonormal eigenfunctions of the undamped wave operator.

6. Contributors

**Lassi Paunonen: **

Project leader, developer (2017–)

**Mikko Aarnos: **

Developer (2018)

**Jukka-Pekka Humaloja: **

Developer (2018–)

Appendix A External Links

• •

https://sysgrouptampere.wordpress.com — Systems Theory Research Group at Tampere University.

• •

https://github.com/lassipau/rorpack/ — RORPack at GitHub

• •

https://lassipau.github.io/rorpack/ — RORPack homepage (hosted by GitHub Pages)

• •

http://mathesaurus.sourceforge.net/matlab-numpy.html — Useful information on differences in the Matlab and Python/Numpy syntax for vectors and matrices.

• •

https://chebfun.org — Homepage of the Chebfun Project.