Model Predictive Control of Parabolic PDE Systems under Chance Constraints

Abstract

Model predictive control (MPC) heavily relies on the accuracy of the system model. Nevertheless, process models naturally contain random parameters. To derive a reliable solution, it is necessary to design a stochastic MPC. This work studies the chance constrained MPC of systems described by parabolic partial differential equations (PDEs) with random parameters. Inequality constraints on time- and space-dependent state variables are defined in terms of chance constraints. Using a discretization scheme, the resulting high-dimensional chance constrained optimization problem is solved by our recently developed inner–outer approximation which renders the problem computationally amenable. The proposed MPC scheme automatically generates probability tubes significantly simplifying the derivation of feasible solutions. We demonstrate the viability and versatility of the approach through a case study of tumor hyperthermia cancer treatment control, where the randomness arises from the thermal conductivity coefficient characterizing heat flux in human tissue.

1. Introduction

Partial differential equations (PDEs) are models commonly used for describing spatial and/or temporal variations in natural, physical, and industrial processes [1,2]. Conventionally, the parameters involved in PDEs, e.g., thermal conductivity, diffusivity, viscosity, charge density, and material density [1,3], are taken to be deterministic coefficients that are sometimes prescribed from experiences or fixed through experimental measurements. Nevertheless, such nominal values do not precisely reflect the nature of the parameters, since they fluctuate in the real process and, hence, are uncertain. In addition, the operational environment of the physical system and uncontrolled external sources, e.g., ambient temperature and pressure, can also impose a significant impact. These uncertainties constitute what is known as input uncertainties (or random data), which propagate through the system and incur discrepancies between the real (actual) responses and the outputs predicted by the nominal PDE model.

Various simulation-based uncertainty quantification (UQ) methods have been developed for capturing propagated uncertainties from the inputs to the outputs (state variables) of a random PDE system [4,5,6,7,8,9,10,11,12,13,14,15,16,17]. Such studies have recently triggered various investigations on the mathematical analysis, numerical solution methods, and applications of PDE optimization problems with input uncertainties [11,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33]. In particular, studies on the state-constrained optimization of PDE systems with random data constitute the state-of-the-art. For such problems, an appropriate definition of inequality constraints of random state variables is basically a perquisite to render the optimization problem as well-defined. Hence, recent studies [19,20,30,33] considered chance constrained PDE optimization problems, i.e., to impose chance constraints on the state variables. However, these studies were focused more on elliptic PDE systems. Owing to the wide-spread use of parabolic PDEs in various engineering applications with distributed parameter systems and time-dependent behaviors, studies on optimization problems with state-constrained parabolic PDE systems involving random data are necessary.

Model predictive control (MPC) provides advanced control strategies for systems with complex state-space models involving control and state constraints. Recent studies show that the MPC scheme has further potential uses in engineering applications with distributed parameter models, such as parabolic PDEs, e.g., thermal processes [34,35,36] and process engineering applications [37,38,39,40,41,42], etc. Nevertheless, these studies investigated deterministic MPC without considering uncertainties. MPC for state-constrained distributed parameter systems with random data has not yet been well-studied.

Since reliable predictions are heavily dependent on the accuracy of the underlying model, a deterministic MPC scheme is bound to lead to infeasible state trajectories with inappropriately designed control strategies. Therefore, this work studies the model predictive control of processes described by parabolic PDE systems with random parameters. We propose a chance constrained MPC scheme for the control of such systems with random coefficients by imposing probabilistic measures for the reliable satisfaction of the state constraints.

The MPC of parabolic PDE systems with chance constrained state variables calls for the solution of a chance constrained PDE optimization problem on each prediction horizon. For a numerical solution, we transform it through an appropriate discretization into a high dimensional chance constrained optimization problem. Although many studies have been made to develop solution approaches [43,44,45,46], open challenges remain in particular for addressing large-scale complex problems. In general, there are three major difficulties in solving chance constrained problems on finite or infinite dimensional spaces: (a) the probability function of a chance constraint can be non-convex and non-differentiable, (b) this function is usually hard to evaluate, and (c) the feasibility to the chance constraint can be hard to verify. Even the knowledge of convexity and differentiability may not ensure a simpler evaluation of chance constraints. Therefore, for complex models, chance constraints can only be evaluated through approximation methods.

In this study, we focus on the development of a numerical solution approach to chance constrained MPC of semi-linear parabolic PDE systems. Hence, to circumvent potential non-smoothness and inherent computational difficulties associated with chance constraints, we extend the applicability of our smoothing inner–outer approximation method [47,48] to chance constrained MPC of such PDE systems. The proposed approach entails two important advantages. First, on each prediction horizon, the chance constrained problem is approximated by two computationally amenable smooth (inner and outer) optimization problems. The respective feasible sets and optimal solutions of the smoothing problems are convergent to the feasible set and solutions of the chance constrained problem (see [20,47,48] for detailed theoretical analysis). Second, the approximation scheme automatically generates probability tubes, whereby the inner approximation guarantees feasibility of the predicted state trajectories to the chance constrained MPC problem. In this way, the effort for constructing a probability tube, which is the case in traditional literature concerning tube-based stochastic MPC [49,50,51,52,53,54], is avoided.

One of the important areas to deal with challenging PDE models is medicine, e.g., physiologically based pharmacokinetic modeling for drug development [55], hyperthermia treatment in cancer therapy [56,57,58] and medical image analysis [59,60]. To demonstrate our approach, we extend the PDE-constrained optimal control problem of hyperthermia treatment by considering uncertain parameters, imposing chance constraints, and applying the receding horizon strategy for computing optimal as well as reliable control strategies for cancer therapy.

2. Problem Statement and Assumptions

2.1. Problem Statement

We consider a chance constrained model predictive controller (CCMPC) of a semi-linear parabolic PDE system:

$$\begin{array}{ccc}& & \left(\mathrm{CCMPC}\right)\hfill \end{array}$$

$$\begin{array}{ccc}& & \underset{u}{min}\left\{J\left(u\right)=E\left[{\int }_{{\widehat{t}}_z}^{{\widehat{t}}_z+H}\left\{\parallel y(u,\xi ;t,·)-y_d{(t,·)\parallel }_{L^2\left(D\right)}^2+\frac{\lambda }{2}{\parallel u(t,·)\parallel }_{L^2\left(D\right)}^2\right\}dt\right]\right\}\hfill \\ & & \mathrm{subject}\mathrm{to}:\hfill \end{array}$$

(1)

$$\begin{array}{ccc}& & \frac{\partial y}{\partial t}-{\nabla }_x\left[\kappa \left(\xi \right){\nabla }_{x}y\right]=f(u,t,x)\mathrm{in}Q\times \mathsf{\Omega },\hfill \end{array}$$

(2)

$$\begin{array}{ccc}& & -{\kappa }_0{\nabla }_{x}y·\mathbf{n}=g(y,t,x)\mathrm{on}({\widehat{t}}_z,{\widehat{t}}_z+H]\times \partial D\times \mathsf{\Omega },\hfill \end{array}$$

(3)

$$\begin{array}{ccc}& & y({\widehat{t}}_z,x)=y_{{\widehat{t}}_z}\left(x\right)\mathrm{in}D,\hfill \end{array}$$

(4)

$$\begin{array}{ccc}& & Pr\left\{y(u,\xi ;t,x)\le y_{max}\right\}\ge \alpha ,\mathrm{in}Q,\hfill \end{array}$$

(5)

$$\begin{array}{ccc}& & u_{min}\le u(t,x)\le u_{max},\mathrm{in}Q,\hfill \end{array}$$

(6)

$$\begin{array}{ccc}& & {\widehat{t}}_{z+1}:={\widehat{t}}_z+\Delta t,z=1,\dots ,N_t^{\prime },\hfill \end{array}$$

(7)

where $D\subset {\mathbb{R}}^n$ ($n=1,2,3$) is a bounded spatial domain with boundary $\partial D$; $x\in D$ represents the spatial position; ∇ denotes derivative with respect to x; $\mathbf{n}$ stands for a unit outward normal vector to $\partial D$; the interval $({\widehat{t}}_z,{\widehat{t}}_z+H]$ is a prediction time-horizon with a finite length H; $Q:=({\widehat{t}}_z,{\widehat{t}}_z+H]\times D$; ${\xi }^{\top }\in \mathsf{\Omega }\subset {\mathbb{R}}^p$ is a vector of random input variables from a complete probability space $\left(\mathsf{\Omega },\Sigma ,Pr\right)$ with $\sigma$-algebra $\Sigma$ and probability measure $Pr(·)$ associated with a given continuous probability density function $\varphi \left(\xi \right)$, which can be Gaussian or non-Gaussian distributed. Moreover, $E[·]$ represents the expected value operator with respect to $\varphi \left(\xi \right)$.

$\kappa \left(\xi \right)$ is considered as a random coefficient in the PDE system with nominal value ${\kappa }_0=E\left[\kappa \left(\xi \right)\right]$, where $\xi$ represents model parameter uncertainty. The system is driven by a distributed forcing term $f(u,t,x)$ with a control variable u. Equation (3) describes the interaction of the system with the surrounding environment through the boundary $\partial D$, which is specified by the boundary function $g(y,t,x)$.

The control variable u is deterministic and belongs to the set $U=\{u\in L^2\left(Q\right)|u_{min}\le u\le u_{max}\}$. The state variable y depends on the control u and the random parameter $\xi$, which is expressed as $y(u,\xi ;t,x)$, $(t,x)\in Q$. Since y is a random variable, the chance constraint (5) specifies the probability of holding restrictions on the state variable, where $\alpha \in \left[0,1\right]$ is a user-specified level of reliability. Since the chance constraint (5) is required to hold point-wise over the spatial domain D at each time instant $t\in ({\widehat{t}}_z,{\widehat{t}}_z+H]$, this formulation has the advantage that the user can selectively specify higher or lower reliability levels at given locations x and instant of time t.

Associated to the CCMPC defined above, we further define

• the probability function

  $$p(u;t,x):=Pr\left\{y(u,\xi ;t,x)\le y_{max}\right\},$$
• feasible set

  $$\mathcal{P}=\{u\in U|p(u;t,x)\ge \alpha ,(t,x)\in Q\}.$$

2.2. Assumptions

Although formal mathematical analytic investigations are not provided here, the discussion needs to be framed in an appropriate function space. Moreover, to render the problem as well-posed as well as to validate the numerical solution approach, the problem setting needs to satisfy standard assumptions and the relevant function spaces need to be defined.

• Let $L^2\left(Q\right)$ be the space of square-integrable functions over Q, then define

  $$\mathcal{F}:=\left\{v:\mathsf{\Omega }\to L^2{\left(Q\right),v\mathrm{is}\mathrm{measurable},\parallel v\parallel }_{\mathcal{F}}:={\int }_{\mathsf{\Omega }}{\parallel v(\xi ,·)\parallel }_{L^2\left(Q\right)}\varphi \left(\xi \right)d\xi <+\infty \right\}.$$
• Let $L^2\left(S\right)$ be the space of square-integrable functions over $S:=({\widehat{t}}_z,{\widehat{t}}_z+H]\times \partial D$, then define

  $$\mathcal{G}:=\left\{v:\mathsf{\Omega }\to L^2{\left(S\right),v\mathrm{is}\mathrm{measurable},\parallel v\parallel }_{\mathcal{G}}:={\int }_{\mathsf{\Omega }}{\parallel v(\xi ,·)\parallel }_{L^2\left(S\right)}\varphi \left(\xi \right)d\xi <+\infty \right\}.$$
• Using the Sobolev space $H^1\left(D\right)$ and the Hilbert space $\mathcal{W}:=L^2\left(\mathsf{\Omega },H^1\left(D\right)\right)$

  $$\mathcal{V}:=\left\{v\in L^2\left(\mathsf{\Omega },H^1\left(D\right)\right)\left|{\parallel v\parallel }_{\mathcal{W}}^2={\int }_{\mathsf{\Omega }}{\int }_D\left(v^2+\kappa {|\nabla v|}^2\right)dx\varphi \left(\xi \right)d\xi <+\infty \right.\right\}.$$

The following assumptions are standard in the literature for PDE systems with a random parameter.

A1: 

Uniform ellipticity condition. There are positive constants ${\kappa }_{max},{\kappa }_{min}$ with ${\kappa }_{max}>{\kappa }_{min}$, such that the random coefficient $\kappa$ satisfies

$${\kappa }_{min}\le \kappa \left(\xi \right)\le {\kappa }_{max},\mathrm{almost}\mathrm{everywhere}.$$

A2: 

For each fixed $u\in U$, the forcing term f has the property that $f(u,·,·)\in L^2\left(Q\right)$.

A3: 

For each fixed $u\in U$, the boundary function $g(u,\xi ;t,x)=g\left(y\right(u,\xi ;t,x);t,x)$ has the property that $g(u,·,·,·)\in \mathcal{G}$.

Here, A1 is a realistic assumption from a practical point of view. It is a crucial assumption also to guarantee that the underlying differential operator of the PDE system is uniformly elliptic. Thereby, under the function spaces defined and for the problem settings specified in assumptions A1–A3, the parabolic PDE system Equation (2) with boundary (3) and initial Equation (4) can be guaranteed to have a unique weak solution, which is continuously dependent on f and g [61].

Hence, we assume that, for each pair $(u,\xi )\in Q\times \mathsf{\Omega }$, the parabolic PDE system Equations (2)–(4) hve a unique solution $y(u,\xi ;t,x),(t,x)\in Q$, so that $y(·,·;t,x),(t,x)\in Q$ is continuous with respect to $(u,\xi )$. Moreover, it is assumed that $y(u,·;t,·)\in L^p(\mathsf{\Omega };\mathcal{V})$, where $\mathcal{V}$ is a Banach space, e.g., $\mathcal{V}=L^2\left(D\right)$ or $\mathcal{V}=H_0^1\left(D\right)$, and

$${\parallel y(u,·,t,·)\parallel }_{L^p(\mathsf{\Omega };\mathcal{V})}^p={\int }_{\mathsf{\Omega }}\left[{\int }_D{\left|y(u,\xi ;t,x)\right|}^{p}dx\right]\varphi \left(\xi \right)d\xi =E\left[{\parallel y(t,·,\xi )\parallel }_{\mathcal{V}}^p\right]<\infty$$

for each $t>0$.

Over each prediction horizon of the MPC scheme, the chance constrained optimal control problem (CCPDE) with the random parabolic PDE system Equations (2)–(4) is to be solved based on the deterministic initial condition Equation (4). The solution yields an optimal control sequence $u^{*}$ over $({\widehat{t}}_z,{\widehat{t}}_z+H]$, and an optimal feedback is obtained by applying the first element of $u^{*}$ to the PDE system. To employ the receding horizon strategy, it is necessary to update the initial time instance ${\widehat{t}}_z$ of the prediction horizon by Equation (7), where $\Delta t$ is a sampling time, and $N_t^{\prime }$ is the number of time discretization points in the operational period $({\widehat{t}}_1,{\widehat{t}}_1+H^{\prime }]$ of the finite length $H^{\prime }$.

3. Solution Approach to CCMPC

3.1. Finite Dimensional Representation

The CCMPC of parabolic PDE systems calls for the solution of a CCPDE optimization problem on each prediction horizon. The theoretical investigation of this problem can be worked out in a more general context of Banach spaces (see [20]).

Moreover, based on the discussion in Section 2, we assume that for each given admissible control u and realization of the random variable $\xi$, the PDE system of Equations (2)–(4) can be solved to obtain $y(u,\xi ;t,x)$. As a result, on each prediction horizon, the problem reduces to

$$\begin{array}{ccc}& & \left({\mathrm{CCPDE}}_r\right)\hfill \end{array}$$

$$\begin{array}{ccc}& & \underset{u}{min}\left\{J\left(u\right)=E\left[{\int }_{{\widehat{t}}_z}^{{\widehat{t}}_z+H}\left\{\parallel y(u,\xi ;t,·)-y_d{(t,·)\parallel }_{L^2\left(D\right)}^2+\frac{\lambda }{2}{\parallel u(t,·)\parallel }_{L^2\left(D\right)}^2\right\}dt\right]\right\}\hfill \\ & & \mathrm{subject}\mathrm{to}:\hfill \end{array}$$

(8)

$$\begin{array}{ccc}& & p(u;t,x):=Pr\left\{y(u,\xi ;t,x)\le y_{max}\right\}\ge \alpha ,(t,x)\in Q,\hfill \end{array}$$

(9)

$$\begin{array}{ccc}& & u\in U.\hfill \end{array}$$

(10)

For numerical computations, the reduced formulation ${\mathrm{CCPDE}}_r$ Equations (8)–(10) need to be transformed into a finite dimensional representation on the prediction horizon through an appropriate space-time discretization technique. Suppose $\{{\widehat{t}}_z,{\widehat{t}}_z+\Delta t,\dots ,{\widehat{t}}_z+(N_t-1)\Delta t\}$ is a set of appropriately chosen time instances, where the number of time discretization points $N_t=[H\Delta t]+1$, and $\{x_1,x_2,\dots ,x_{N_x}\}$ are spatial discretization points generated according to a given finite element method from the domain D, respectively. Subsequently, these yield a vector of discretized controls $\mathbf{u}$ Equation (11), a corresponding admissible set ${\mathcal{U}}_d$ Equation (12), a discretized cost functional $J\left(\mathbf{u}\right)$ Equation (13), and a set of chance constraints $p_{j,k}$ Equation (14) corresponding to each time instance and space point, as follows

$$\begin{array}{ccc}& & \mathbf{u}={[u(t_1,x_1),u(t_2,x_2),\dots ,u(t_{N_t},x_{N_x})]}^T\in {\mathbb{R}}^{N_t\times N_x},\hfill \end{array}$$

(11)

$$\begin{array}{ccc}& & {\mathcal{U}}_d=\{\mathbf{u}\in {\mathbb{R}}^{N_t\times N_x}|\underline{\mathbf{u}}\le \mathbf{u}\le \overline{\mathbf{u}}\},\hfill \end{array}$$

(12)

$$\begin{array}{ccc}& & J\left(\mathbf{u}\right)=E\left[\sum _{j=1}^{N_t}\sum _{k=1}^{N_x}\left\{{\left(y(u(t_j,x_k),\xi ,t_j,x_k)-y_d(t_j,x_k)\right)}^2+\frac{\lambda }{2}u{(t_j,x_k)}^2\right\}\right],\hfill \end{array}$$

(13)

$$\begin{array}{ccc}& & p_{j,k}\left(\mathbf{u}\right)=Pr\{y(\mathbf{u},\xi ,t_j,x_k)-y_{max}(t_j,x_k)\le 0\}\ge \alpha ,j=1,\dots ,N_t;k=1,\dots ,N_x,\hfill \end{array}$$

(14)

where $\underline{\mathbf{u}}={[u_{min},\dots ,u_{min}]}^{\top }\in {\mathbb{R}}^{N_t\times N_x}$, $\overline{\mathbf{u}}={[u_{max},\dots ,u_{max}]}^{\top }\in {\mathbb{R}}^{N_t\times N_x}$. Consequently, a finite dimensional representation of CCPDE${}_r$ takes the form

$$\begin{array}{ccc}\hfill \left(CCOPT\right)& & \underset{\mathbf{u}\in {\mathcal{U}}_d}{min}J\left(\mathbf{u}\right)\hfill \\ & & \mathrm{subject}\mathrm{to}:\hfill \end{array}$$

(15)

$$\begin{array}{ccc}& & p_{j,k}\left(\mathbf{u}\right)=Pr\{s_{j,k}(\mathbf{u},\xi )\le 0\}\ge \alpha ,j=1,\dots ,N_t;k=1,\dots ,N_x,\hfill \end{array}$$

(16)

where $s_{j,k}(\mathbf{u},\xi ):=y(\mathbf{u},\xi ,t_j,x_k)-y_{max}(t_j,x_k)$ and with the feasible set

$$\mathcal{P}=\left\{u\in {\mathcal{U}}_d|p_{j,k}\left(\mathbf{u}\right)\ge \alpha ,j=1,\dots ,N_t;k=1,\dots ,N_x\right\}.$$

(17)

As mentioned above, the major computational hurdle associated with chance constrained optimization problems is the difficulty in evaluating the probability

$$p_{j,k}\left(u\right)=Pr\left\{s_{j,k}(\mathbf{u},\xi )\le 0\right\}=\underset{\left\{\xi \in \mathsf{\Omega }|s_{j,k}(\mathbf{u},\xi )\le 0\right\}}{\int }\varphi \left(\xi \right)d\xi .$$

(18)

The probability function $p_{j,k}\left(u\right)$ usually does not have a closed form analytic representation. Except for special types of sets $\left\{\xi \in \mathsf{\Omega }|s_{j,k}(\mathbf{u},\xi )\le 0\right\}$, it is generally difficult to directly evaluate the integral on the right hand-side of Equation (18) to obtain the value of $p_{j,k}\left(u\right)$ for a given u. Furthermore, $p_{j,k}\left(u\right)$ can be non-differentiable, non-convex, and commonly guaranteeing the feasibility of an optimal solution is not trivial. Moreover, the theoretical assumptions for the differentiability of $p_{j,k}\left(u\right)$ are quite restrictive and also the knowledge of its differentiability furnishes little computational advantages. Therefore, to overcome these difficulties, we use our recently proposed smoothing inner–outer approximation strategy [47]. The viability of this approximation strategy has been demonstrated for both finite [47,48] and infinite dimensional [20] CCOPT problems. Here, we extend the inner–outer approach for the solution of chance constrained MPC problems for parabolic PDE systems. Using this method, we guarantee the feasibility of the state trajectories in MPC, despite the presence of uncertainties and, in particular, without requiring an extra effort for constructing a probability tube.

3.2. Inner-Outer Approximation Strategy

Here, we briefly summarize and adapt the inner–outer approximation approach of Geletu et al. [47] to a CCMPC problem, highlighting its relevant properties.

3.2.1. The Smoothing Approximation Problems

By dropping the indices from Equation (16) for the sake of brevity, consider a single chance constraint

$$p\left(u\right)=Pr\left\{s\right(u,\xi )\le 0\}\ge \alpha .$$

Define two parametric functions $\psi (\tau ,u)$ and $\phi (\tau ,u)$, with $\tau \in (0,1)$, based on the Geletu–Hoffmann function $\mathsf{\Theta }(\tau ,s)$ (Geletu et al. [47], also [20,48]) such that

$$\psi (\tau ,u)=E\left[\mathsf{\Theta }(\tau ,s(u,\xi \left)\right)\right]\mathrm{and}\phi (\tau ,u)=E\left[\mathsf{\Theta }(\tau ,-s(u,\xi \left)\right)\right],$$

where

$$\mathsf{\Theta }(\tau ,s)=\frac{1+m_1\tau }{1+m_2\tau exp\left(-\frac{s}{\tau }\right)},$$

(19)

and $m_1,m_2$ are constants with $0<m_2\le m_1/(1+m_1)$.

Now, we define two parametric optimization problems in

Table 1. Inner and outer approximation problems.

Inner-approximation problem	Outer-approximation problem
$\begin{array}{cc}\hfill {\left(IA\right)}_{\tau }& \underset{u}{min}J\left(u\right)\hfill \\ \hfill & \mathrm{subject}\mathrm{to}:\hfill \\ \hfill & \psi (u,\tau )\le 1-\alpha ,\hfill \\ & u\in {\mathcal{U}}_d,\hfill \end{array}$	$\begin{array}{cc}\hfill {\left(OA\right)}_{\tau }& \underset{u}{min}J\left(u\right)\hfill \\ \hfill & \mathrm{subject}\mathrm{to}:\hfill \\ \hfill & \phi (u,\tau )\ge \alpha ,\hfill \\ & u\in {\mathcal{U}}_d,\hfill \end{array}$

with

$$M\left(\tau \right)=\{u\in {\mathcal{U}}_d|\psi (u,\tau )\le 1-\alpha \}\mathrm{and}S\left(\tau \right)=\{u\in {\mathcal{U}}_d|\phi (u,\tau )\ge \alpha \}$$

are feasible sets of the problems (IA)${}_{\tau }$ and (OA)${}_{\tau }$, respectively.

3.2.2. Properties of the Approximations

• Properties of the Approximation Functions

  (a)

  Uniform limit. Property 2 of Theorem 3.6 in [47] implies that both functions $1-\psi (\tau ,u)$ and $\phi (\tau ,u)$ closely resemble the function $p\left(u\right)$, as the parameter $\tau$ moves closer to zero from above. This uniform approximation property has two important consequences.

  (${\mathrm{a}}_1$)

  $p(·)$ is a continuous function (Corollary 3.8 in [20]).

  (${\mathrm{a}}_2$)

  The function $\psi (\tau ,·)$ and $\phi (\tau ,·)$ are differentiable with regard to u (Theorem 3.6 in [47]).

  Note that ${\mathrm{a}}_1$ and ${\mathrm{a}}_2$ hold true irrespective of the differentiability of the probability function $p(·)$. Therefore, a gradient-based optimization algorithm can be used to solve the smooth approximation problems $\mathrm{I}{\mathrm{A}}_{\tau }$ and $\mathrm{O}{\mathrm{A}}_{\tau }$, respectively, to obtain approximate optimal solutions to the (generally non-smooth) CCOPT problem.

  (b)

  If the same assumptions as in [62] hold true, then the probability function $p(·)$ is differentiable and the gradients $\nabla (1-\psi (\tau ,u\left)\right)$ and $\nabla \phi (\tau ,u)$ can be shown to converge, with regard to $\tau$, to the gradient $\nabla p\left(u\right)$ (Theorem 16 in [47] and Theorem 4.10 in [48]).

• Properties of the inner–outer approximation problems
  Recall that the feasible set of CCOPT is given by $\mathcal{P}=Pr\left\{s\right(u,\xi )\le 0\}\ge \alpha .$

  (a)

  Inner–outer approximation. The feasible set $M\left(\tau \right)$ is always a subset (inner approximation) of the feasible set $\mathcal{P}$ of CCOPT, while the feasible set $S\left(\tau \right)$ is a superset (outer approximation) of the feasible set $\mathcal{P}$ of CCOPT; i.e.,

  $$\mathcal{M}\left(\tau \right)\subset \mathcal{P}\subset S\left(\tau \right),\mathrm{for}\mathrm{any}\tau \in (0,1).$$

  (b)

  Monotonicity of the inner–outer approximations to ensure the enclosure of the feasible set:

  $$\mathcal{M}\left({\tau }_2\right)\subset \mathcal{M}\left({\tau }_1\right)\subset \mathcal{P}\subset S\left({\tau }_1\right)\subset S\left({\tau }_2\right),\mathrm{for}0<{\tau }_1<{\tau }_2<1.$$

  (c)

  Convergence of the feasible sets of the approximations:

  $$\underset{\tau \searrow 0^{+}}{lim}\mathcal{M}\left(\tau \right)=\mathcal{P},\underset{\tau \searrow 0^{+}}{lim}S\left(\tau \right)=\mathcal{P},$$

  where the first convergence needs the regularity condition $cl\left\{u\in U\left|p\right(u)>\alpha \right\}=\mathcal{P}$ (see Geletu et al. [47]).

  (d)

  Obviously, $J_{IA}(·)$ of ${\left(IA\right)}_{\tau }$ is an upper bound to the objective function of ${\left(OA\right)}_{\tau }$ and converge to each other as $\tau$ moves closer to 0 from above.

• The algorithm for computation
  The following algorithm presents the computation scheme of the inner–outer approximation approach for solving a CCOPT problem.

In Algorithm 1, $J_{IA}\left(u_{{\tau }_k}^{*}\right)$ and $J_{OA}\left({\tilde{u}}_{{\tau }_k}^{*}\right)$ represent the optimal objective function values of ${\left(\mathrm{I}\mathrm{A}\right)}_{{\tau }_k}$ and ${\left(\mathrm{O}\mathrm{A}\right)}_{{\tau }_k}$ at the respective optimal points $u_{{\tau }_k}^{*}$ and ${\tilde{u}}_{{\tau }_k}^{*}$. The algorithm stops when the difference of optimal values of the inner and outer approximations is less than a given threshold.

Algorithm 1: Computation scheme for the inner–outer approximation method.

1: Choose an initial parameter ${\tau }_0\in (0,1)$; 2: Solve the optimization problems ${\left(\mathrm{IA}\right)}_{{\tau }_0}$ and ${\left(\mathrm{OA}\right)}_{{\tau }_0}$; 3: Select the termination tolerance tol; 4: Set k←0; 5: while ($|J_{IA}\left(u_{{\tau }_k}^{*}\right)-J_{OA}\left({\tilde{u}}_{{\tau }_k}^{*}\right)|$ > tol) do 6: Reduce the parameter τk (e.g., τk+1 = ρτk, for ρ ∈ (0, 1)); 7:  Set k ← k + 1; 8: Solve the optimization problems ${\left(\mathrm{IA}\right)}_{{\tau }_k}$ and ${\left(\mathrm{OA}\right)}_{{\tau }_k}$; 9: end while

Remark 1.

Referring the chance constraint $Pr\left\{s(u,\xi ;t,x)\le 0\right\}\ge \alpha ,inQ=(\widehat{t},\widehat{t}+H]\times D$ of CCMPC with $s(u,\xi ;t,x)=y(u,\xi ;t,x)-y_{max}$, define

$$\psi (\tau ,u;t,x)=E\left[\mathsf{\Theta }(\tau ,s(u,\xi ;t,x\left)\right)\right]and\phi (\tau ,u;t,x)=E\left[\mathsf{\Theta }(\tau ,-s(u,\xi ;t,x\left)\right)\right].$$

(20)

Using the same parametric function $\mathsf{\Theta }(\tau ,·)$ as in Equation (19) and let $\mathcal{P}\left(t\right)=\{u\in {\mathcal{U}}_d|p(u;t,x)\ge \alpha ,x\in D\}$,$M(\tau ;t)=\{u\in {\mathcal{U}}_d|\psi (u,\tau ;t,x)\le 1-\alpha ,x\in D\}$ and $S(\tau ;t)=\{u\in {\mathcal{U}}_d|\phi (u,\tau ;t,x)\ge \alpha ,x\in D\}$,$t\in (\widehat{t},\widehat{t}+H]$. It follows that,

$$M(\tau ;t)\subset \mathcal{P}\left(t\right)\subset S(\tau ;t),t\in (\widehat{t},\widehat{t}+H].$$

This implies that the feasible sets of the inner and outer approximations automatically define deterministic inner and outer enveloping tubes to the feasible set of the chance constraint over each prediction horizon  $(\widehat{t},\widehat{t}+H]$ (see Figure 1). These tubes guarantee the feasibility of the solution on each prediction horizon and will be made tighter by decreasing parameter τ by $\Delta \tau$, from one prediction horizon to the next, yielding robust constraint satisfaction. The tubes are constructed automatically, which neither incurs additional computation expenses nor does it require additional assumptions to simplify the problem.

4. Numerical Implementation

The CCMPC of the parabolic PDE system is reformulated as CCPDE with time discretization, where on each prediction horizon a CCPDE problem is solved by updating the initial condition Equation (4) at the time moment ${\widehat{t}}_z$ following the receding horizon strategy Equation (7). For this purpose, we use the backward difference for the time discretization with $N_t$ time instances and finite element space discretization of the PDE system with $N_x$ space points, respectively. N quasi-Monte Carlo samples are generated for $\xi$. As a result, we obtain two finite dimensional NLP problems for the inner and outer approximations as follows:

$$\begin{array}{ccc}\hfill {\left({\mathrm{NLP}}_{IA}\right)}_{\tau }& & \underset{u}{min}\left\{\frac{1}{N}\sum _{i=1}^N\sum _{j=1}^{N_t}\sum _{k=1}^{N_x}\left[{\left(y(u(t_j,x_k),{\xi }_i,t_j,x_k)-y_d\right)}^2+\frac{\lambda }{2}u{(t_j,x_k)}^2\right]\right\}\hfill \\ & & \mathrm{subject}\mathrm{to}:\hfill \end{array}$$

(21)

$$\begin{array}{ccc}& & \mathsf{\Psi }(\tau ,u(t_j,x_k))\le 1-\alpha ,j=1,\dots ,N_t;k=1,\dots ,N_x,\hfill \end{array}$$

(22)

$$\begin{array}{ccc}& & u_{min}\le u(t_j,x_k)\le u_{max},j=1,\dots ,N_t;k=1,\dots ,N_x,\hfill \end{array}$$

(23)

$$\begin{array}{ccc}\hfill \mathrm{where}& & \mathsf{\Psi }(\tau ,u(t_j,x_k))=\frac{1}{N}\sum _{i=1}^N\frac{1+m_1\tau }{1+m_2\tau exp(y_{max}-y(u(t_j,x_k),{\xi }_i,t_j,x_k))},\hfill \end{array}$$

(24)

$$\begin{array}{ccc}\hfill {\left({\mathrm{NLP}}_{OA}\right)}_{\tau }& & \underset{u}{min}\left\{\frac{1}{N}\sum _{i=1}^N\sum _{j=1}^{N_t}\sum _{k=1}^{N_x}\left[{\left(y(u(t_j,x_k),{\xi }_i,t_j,x_k)-y_d\right)}^2+\frac{\lambda }{2}u{(t_j,x_k)}^2\right]\right\}\hfill \\ & & \mathrm{subject}\mathrm{to}:\hfill \end{array}$$

(25)

$$\begin{array}{ccc}& & \mathsf{\Phi }(\tau ,u(t_j,x_k))\ge \alpha ,j=1,\dots ,N_t;k=1,\dots ,N_x,\hfill \end{array}$$

(26)

$$\begin{array}{ccc}& & u_{min}\le u(t_j,x_k)\le u_{max},j=1,\dots ,N_t;k=1,\dots ,N_x,\hfill \end{array}$$

(27)

$$\begin{array}{ccc}\hfill \mathrm{where}& & \mathsf{\Phi }(\tau ,u(t_j,x_k))=\frac{1}{N}\sum _{i=1}^N\frac{1+m_1\tau }{1+m_2\tau exp(y(u(t_j,x_k),{\xi }_i,t_j,x_k)-y_{max})}.\hfill \end{array}$$

(28)

Numerically, the spatial discretization of the PDEs is realized by the finite element method in the FEniCS library (FEniCS project [63,64]). Alternatively, libraries supporting spectral discretization could be used. Spectral discretization can improve the numerical accuracy of discretization, allowing larger intervals and leading to higher computational efficiency. However, in comparison to such libraries, FEniCS also allows the solution of the PDE-constrained optimization problems through the Dolfin-Adjoint package [65]. Despite the automated finite element spatial discretization provided by FEniCS, Dolfin-Adjoint does not directly support the solution of state-constrained optimization problems with PDEs. Therefore, to solve the problems NLP${}_{IA}$ and NLP${}_{OA}$ under FEniCS, we embed constraints (22) and (26) into the respective objective functions. In particular, when state variables are involved in inequality constraints explicitly or implicitly, the Moreau–Yosida regularization [66,67] can render the problem amenable for solving using FEniCS. Thus, we use the smoothed version of the Moreau–Yosida regularization and express the regularized problems as

$$\begin{array}{ccc}& & {\left({\mathrm{RNLP}}_{IA}\right)}_{\tau }\hfill \end{array}$$

$$\begin{array}{ccc}& & \underset{u}{min}\left\{\frac{1}{N}\sum _{i=1}^N\sum _{j=1}^{N_t}\sum _{k=1}^{N_x}\left[{\left(y(u(t_j,x_k),{\xi }_i,t_j,x_k)-y_d\right)}^2+\frac{\lambda }{2}u{(t_j,x_k)}^2+\frac{\gamma }{2}{R}_{IA}{(\tau ,u(t_j,x_k))}^2\right]\right\}\hfill \\ & & \mathrm{subject}\mathrm{to}:\hfill \end{array}$$

(29)

$$\begin{array}{ccc}& & u_{min}\le u(t_j,x_k)\le u_{max},j=1,\dots ,N_t;k=1,\dots ,N_x,\hfill \end{array}$$

(30)

$$\begin{array}{ccc}& & {\left({\mathrm{RNLP}}_{OA}\right)}_{\tau }\hfill \end{array}$$

$$\begin{array}{ccc}& & \underset{u}{min}\left\{\frac{1}{N}\sum _{i=1}^N\sum _{j=1}^{N_t}\sum _{k=1}^{N_x}\left[{\left(y(u(t_j,x_k),{\xi }_i,t_j,x_k)-y_d\right)}^2+\frac{\lambda }{2}u{(t_j,x_k)}^2+\frac{\gamma }{2}{R}_{OA}{(\tau ,u(t_j,x_k))}^2\right]\right\}\hfill \\ & & \mathrm{subject}\mathrm{to}:\hfill \end{array}$$

(31)

$$\begin{array}{ccc}& & u_{min}\le u(t_j,x_k)\le u_{max},j=1,\dots ,N_t;k=1,\dots ,N_x,\hfill \end{array}$$

(32)

where $\gamma$ is a penalization parameter, and $R_{IA}$ and $R_{OA}$ are (discretized) Moreau–Yosida regularization terms expressed as

$$\begin{array}{ccc}\hfill R_{IA}(\tau ,u(t_j,x_k))& =& \frac{\mathsf{\Psi }(\tau ,u(t_j,x_k))-1+\alpha +\sqrt{{(\mathsf{\Psi }(\tau ,u(t_j,x_k))-1+\alpha )}^2+\chi }}{2},\hfill \end{array}$$

(33)

$$\begin{array}{ccc}\hfill R_{OA}(\tau ,u(t_j,x_k))& =& \frac{\alpha -\mathsf{\Phi }(\tau ,u(t_j,x_k))+\sqrt{{(\alpha -\mathsf{\Phi }(\tau ,u(t_j,x_k)))}^2+\chi }}{2}.\hfill \end{array}$$

(34)

In addition, the Moreau–Yosida regularization strategy inserts the $L_1$ penalty function $max\{0,f(a\left)\right\}$ to the objective function for the constraints of the form $f\left(a\right)\le 0$. However, to avoid working with this non-smooth expression, we use the equality $max\{0,f(a\left)\right\}=\left(f\right(a)+|f\left(a\right)\left|\right)2$ and replace the absolute value by the smoothing approximation $\left|f\left(a\right)\right|\approx \sqrt{f{\left(a\right)}^2+\chi }$ as in Equations (33) and (34), where $\chi \in {\mathbb{R}}^{+}$ is a sufficiently small smoothing parameter.

Figure 2 shows our implementation details. On each prediction horizon, the time and space discretization as well as random sampling are applied to obtain the finite dimensional representation of smoothing inner and outer approximations of the PDE-constrained optimization problem, (RNLP${}_{IA}$)${}_{\tau }$ and (RNLP${}_{OA}$)${}_{\tau }$, respectively. Note that the samples are generated only once and then used to reformulate the initial problem on each prediction horizon. Subsequently, the Moreau–Yosida regularization is employed to regularize (RNLP${}_{IA}$)${}_{\tau }$ and (RNLP${}_{OA}$)${}_{\tau }$. The regularized problems are then solved using Algorithm 1 for decreasing values of the approximation parameter $\tau$ coupled with a sampling strategy for the random variable $\xi$ to evaluate the expected values in Equation (20) and their gradients $\nabla (1-\psi (\tau ,u\left)\right)$ and $\nabla \phi (\tau ,u)$.

We implement ${\mathrm{RNLP}}_{IA}$ and ${\mathrm{RNLP}}_{OA}$ directly in the finite element code and solve them using Dolfin-Adjoint, which provides automatic differentiation and derives adjoint and tangent linear models from the FEniCS-based forward model for the gradient calculation. At each iteration of the optimization algorithm, $y(u,{\xi }_i,x,t_j)$ of the PDE system Equations (2)–(4) will be obtained for each time instant of the random variable ${\xi }_i$. N-times of solutions of the PDE are necessary, since $y(u,{\xi }_i,x,t_j)$ is required to calculate the values and gradients of the cost functional in ${\mathrm{RNLP}}_{IA}$ and ${\mathrm{RNLP}}_{OA}$.

Furthermore, problems ${\mathrm{RNLP}}_{IA}$ and ${\mathrm{RNLP}}_{OA}$ are solved in a double loop: in the first loop over the initial time instance ${\widehat{t}}_z$, updated by Equation (7), and in the second loop over decreasing values of $\tau$. To speed up the computation, we use parallel computing for the finite element code by spreading the solution calculations of the PDE system for each realization of $\kappa \left({\xi }_i\right)$ and by allowing function evaluation and gradient computation over several processor threads through the MPI standard. It is also worth noting that MPI is the only method that allows running such Dolfin-Adjoint finite element code in parallel. Finally, on each prediction horizon, the resulting NLP problems are solved in parallel by the primal-dual interior-point method (the IPOPT library [68]) with the parallel linear solver HSL_MA97 to obtain optimal control sequences. The solution is used as an initial guess for the optimization on the next prediction horizon, which means a so-called warm MPC strategy.

5. Case Study

To demonstrate the proposed approach, we consider here a hyperthermia cancer treatment problem as a case study. Hyperthermia treatment (HT) consists of the heating of malignant tumors to subdue or eradicate the growth of tumor cells from a given organ [56,57]. Depending on the tissue area to be covered, the therapy is divided into three categories: local, regional, and whole-body HT. While the latter two methods act on the whole body or limbs, a local hyperthermia affects only cancerous tissues [69], for instance, a small skin area in case of skin cancer. This allows us to accurately model the tumor domain, formulate an optimization problem, and find the optimal strategy for thermal treatment, which is considered in this case study.

5.1. Mathematical Model

The Pennes bioheat equation [70]

$$\rho ·c_t·\frac{\partial T}{\partial t}-{\nabla }_x\left[\kappa {\nabla }_{x}T\right]+\omega ·c_b·\left(T-T_a\right)=Q\left(u\right),\mathrm{in}(t,t+H]\times D$$

(35)

is a parabolic PDE frequently used to compute the temperature distribution T in the therapy tissue domain D in response to the heating power distribution Q, which, in terms of local HT, can be waves of different frequency ranges and magnetic fields. Equation (35) has the following coefficients and notations: $\rho$—tissue density, $c_t$—specific heat capacity of tissue, $T_a$—arterial blood temperature, $c_b$—specific heat capacity of blood, $\kappa$—tissue thermal conductivity, and $\omega$—perfusion of blood.

5.2. Source of Uncertainty

Equation (35) is widely used in modeling of hyperthermia treatment and often considered to be deterministic. However, its parameters (35) are found to differ from patient to patient as human tissue characteristics are not identical among all individuals [71,72,73]. In addition, laboratory experiments and measurements for a single patient introduce observational error, for instance, owing to the complexity of the measured tissue or organ. Therefore, the model parameters $\rho$; $c_t$; $\kappa$; $\omega$ and $c_b$ are known to be uncertain [73,74]. As a result, these uncertainties propagate through the model to the predicted values of temperature distribution T in the tumor domain and to the amount of heat transferred to healthy neighboring tissues, making them also uncertain [75,76,77] and sensitive to the parameter variations [78].

In this study, we assume the thermal conductivity $\kappa$ to be a time- and space-independent random parameter and its distribution can be obtained from clinical data. This consideration is justified by the fact that in local hyperthermia a relatively small area of tissue is subjected to therapy. As a consequence, the tumor may have properties near to be homogeneous across the domain; i.e., a variation of $\kappa$ in space and time may be insignificant.

5.3. Optimization Problem

Our aim is to regulate the thermal dose through a chance constrained MPC scheme and to keep the tumor tissue $D_t$ uniformly heated to a desired temperature $T_d$ during the treatment (see Figure 3).

At the same time, the heat injection should satisfy the required temperature restriction with a high reliability by avoiding over-heating of the tumor tissue as well as the healthy neighboring tissues $D\backslash D_t$ under all boundary conditions. In this regard, we formulate the following chance constrained HT MPC problem:

$$\begin{array}{ccc}& & \left(\mathrm{CCHTMPC}\right)\hfill \end{array}$$

$$\begin{array}{ccc}& & \underset{u}{min}\left\{E\left[{\int }_{{\widehat{t}}_z}^{{\widehat{t}}_z+H}\left\{\parallel T(u,\xi ;t,·)-T_d{(t,·)\parallel }_{L^2\left(D\right)}^2+\frac{\lambda }{2}{\parallel u(t,·)\parallel }_{L^2\left(D\right)}^2\right\}dt\right]\right\}\hfill \\ & & \mathrm{subject}\mathrm{to}:\hfill \end{array}$$

(36)

$$\begin{array}{ccc}& & \rho ·c_t·\frac{\partial T}{\partial t}-{\nabla }_x\left[\kappa \left(\xi \right){\nabla }_{x}T\right]+\omega ·c_b·\left(T-T_a\right)=W\left(u\right),\mathrm{in}Q\times \mathsf{\Omega },\hfill \end{array}$$

(37)

$$\begin{array}{ccc}& & -{\kappa }_0{\nabla }_{x}T·\mathbf{n}=h\left(T-T_a\right)\mathrm{on}({\widehat{t}}_z,{\widehat{t}}_z+H]\times \partial D\times \mathsf{\Omega },\hfill \end{array}$$

(38)

$$\begin{array}{ccc}& & T({\widehat{t}}_z,x)=T_{{\widehat{t}}_z}\left(x\right)\mathrm{in}D,\hfill \end{array}$$

(39)

$$\begin{array}{ccc}& & Pr\{T(u,\xi ;t,x)\le T_{max}\}\ge \alpha \mathrm{in}Q,\hfill \end{array}$$

(40)

$$\begin{array}{ccc}& & u_{min}\le u(t,x)\le u_{max}\mathrm{in}Q,\hfill \end{array}$$

(41)

$$\begin{array}{ccc}& & {\widehat{t}}_{z+1}:={\widehat{t}}_z+\Delta t,z=1,\dots ,N_t^{\prime },\hfill \end{array}$$

(42)

where the control vector u is bounded by $u_{min}$ and $u_{max}$ and denotes the heat supply to the tumor region $D_t$, the intensity of u is specified by a regularization parameter $\lambda$ > 0. The heat leaving the controlled domain D through the boundary $\partial D$ due to interaction with the environment is determined by the Robin boundary condition Equation (38). h is the heat transfer coefficient, and $T_{{\widehat{t}}_z}$ is the initial temperature on D that sets the initial condition Equation (39). At the initial moment ${\widehat{t}}_1=0$ initial temperature $T(0,x)=T_s$ with $T_s$ as the skin temperature. $T_{max}$ is the maximum allowed temperature distributed on D to avoid overheating of the body, which will be ensured by the upper chance constraint (40).

In previous deterministic studies on HT planning problems, several authors used state constraint $T\ge T_{min}\mathrm{on}D$ (e.g., [76,79,80]), where $T_{min}$ is the minimum allowed therapy temperature. In a stochastic case, it would take the form of lower chance constraint $Pr\{T\ge T_{min}\}\ge \alpha \mathrm{in}Q$. As it will be shown in the next section, this constraint is redundant, since even without it, the temperature holds sufficiently near to $T_d$.

5.4. Parameters and Coefficients

In Equation (37), we set $\rho =c_t=\omega =c_b=1$ for simplicity. The heat transfer coefficient is set to $h={10}^{-5}$. ${\kappa }_0={10}^{-2}$ is used as the nominal value of the thermal conductivity to generate its random values by $\kappa \left({\xi }_i\right)={\kappa }_0{\xi }_i$. The random parameter $\xi$ is described by the four-parameter beta distribution, the probability density function of which ${\beta }_s(x;s,l,a,b)$ has the following relation with the standard two-parameter beta distribution $\beta (x;a,b)$

$${\beta }_s(x;s,l,a,b)=\frac{1}{s}\beta (\frac{x-l}{s},a,b),$$

where s and l stand for the scale and shift parameters, respectively.

Using the quasi-Monte-Carlo method and ${\beta }_s(x;0.2,0.9,2.0,2.0)$, a sample array of $N=300$ samples is generated. Furthermore, to obtain early prediction of constraint violations, to produce corresponding corrective actions, and to ensure numerical stability of the controller, a finite prediction horizon of the length $H=2$ (scaled time) with sampling time $\Delta t=0.025$ is set. This corresponds to $N_t=81$ time instances on each prediction horizon. Usually, in MPC implementation, the computation time for the update should be smaller than the sampling time. However, our numerical study is focused on the investigation of the applicability of the proposed method. In the domain $D=[0,1]\times [0,1]$, we specify $N_x=21\times 21=441$ points and 800 Lagrange elements of the first degree, $40\%$ of which cover the tumor sub-domain $D_t=[0.3,0.7]\times [0.3,0.7]$. As the control is applied only to the sub-domain of the CCHTMPC problem, the final representation has $9\times 9\times 80=6480$ decision variables on each prediction horizon that need to be obtained while the prediction horizon is moving along the operational period (see Figure 4).

The tumor temperature during the HT should be maintained in the range of $[42,44]{}^{\circ }$C. For the numerical calculation, the temperature values in degrees Celsius are scaled to dimensionless quantities. For instance, the maximum allowed temperature of $44{}^{\circ }$C is scaled to $T_{max}=2.05$, the initial temperature, assuming it to be the skin temperature of $32{}^{\circ }$C, is set to $T_s=0$. The desired temperature of $43.7{}^{\circ }$C is linearly interpolated between $T_{max}$ and $T_s$ and corresponds to $T_d=2.0$. Instead of the arterial blood, the capillary blood is considered and assumed to have the temperature of $34.1{}^{\circ }$C, which is scaled to $T_a=0.36$. Other parameters are listed in

Table 2. Optimization parameters.

Parameter		Value
$u_{max}$	Maximum admissible control	30
$u_{min}$	Minimum admissible control	0
$\lambda$	Control regularization parameter	${10}^{-5}$
$\chi$	Moreau–Yosida regularization smoothing parameter	$1.6\times {10}^{-3}$
$\gamma$	Moreau–Yosida regularization penalization parameter	$5.0\times {10}^5$
$\alpha$	Reliability level	0.95
$m_1$, $m_2$	Geletu–Hoffmann function constants	1.0, 0.4

.

6. Results and Discussions

This section presents the results of the numerical experiments. All computations are carried out on a 64-bit Ubuntu Linux machine of version 20.04.2 LTS, equipped with 16 Threads 3.60 GHz Intel Xeon Gold 6244 CPU and 60 GB RAM.

6.1. Results of Optimal Control

For an illustration of the properties of the proposed approach, we first demonstrate the solutions of the inner and outer NLP problems of CCHTMPC only on the first prediction horizon $(0,H]$, which corresponds to the solution of an optimal control problem—CCPDE (i.e., open-loop control without feedback). The optimal control sequences $u_{(0,H]}^{*}(t,x)$ obtained by a decreasing set of the approximation parameter $\tau$ are then applied to the PDE system Equations (37)–(39), leading to optimal temperature distributions $T^{*}(u_{(0,H]}^{*}(t,x),{\xi }_i,t,x),i=1,\dots ,N$; i.e., for the N values of $\kappa \left({\xi }_i\right)$, from which the mean value $E\left[T^{*}(u_{(0,H]}^{*}(t,x),\xi ,t,x)\right]$ is calculated and presented over time in Figure 5, Figure 6 and Figure 7.

The temperature distributions presented in Figure 5 show that the outer problem with a high value of $\tau =0.1$ provides an infeasible solution that exceeds the maximum allowed temperature $T_{max}=2.05$, while the solution of the inner problem remains feasible, but far from the desired temperature $T_d=2.0$. With the decrease of $\tau$, as shown in Figure 5, Figure 6 and Figure 7 and Figure A1 and Figure A2, the solutions of the inner and the outer problem converge to the desired temperature, and with $\tau =5\times {10}^{-5}$ the solution will not violate the chance constraint any more (Figure 7).

The convergence of the inner and outer solutions can be seen more clearly in Figure 8, where the maximum values of mean temperature realizations $E\left[T^{*}(u_{(0,H]}^{*}(t,x),\xi ,t,x)\right]$ are presented along the time and for each value of $\tau$. This Figure 8, similar to Figure 1, shows how the solutions of the outer problem (red curves) narrow to the border of the feasible set of the original problem at $T_{max}=2.05$ due to the constraining effect of the outer deterministic tubes, whereas the solutions of the inner problem (blue curves) expand to the feasible set of the original problem due to the relaxing effect of the inner deterministic tube, as the parameter $\tau$ decreases. It can also be seen that a strong initial heat impulse makes the state constraint prone to violation in the early moments of control $t=[0.0,0.25]$, which could be worsened if a lower reliability level is imposed on the chance constraints. This can be explained by the fact that the minimum allowed temperature $T_{min}=42{}^{\circ }$C is much closer to the maximum allowed temperature $T_{max}=44{}^{\circ }$C than to the initial temperature $T_s=32{}^{\circ }$C, which would force the controller to give a greater thermal momentum to the tumor body at the beginning in order to approach the feasible region $T\ge T_{min}$ as fast as possible.

At the solution, the objective function values obtained under the optimal control for decreasing $\tau$ values are presented in Figure 9. It can be seen, for $\tau =0.01$, the objective function values move close to each other, and for smaller $\tau$ their difference is quite small. This indicates the convergence of the inner and outer solutions. However, $5\times {10}^{-5}$ is the smallest possible value of $\tau$ in the current settings, since there will be a problem of vanishing gradients attributed to the form of the approximation function in Equation (19) when $\tau$ is smaller than this value. This is due to the fact that it moves closer to the Heaviside function for small values of $\tau$ [48]. It should also be noted that the Moreau–Yosida regularization terms Equations (33) and (34) are excluded from the calculation of the optimal objective function values presented in Figure 9, since they contain the parameter $\tau$ and affect the objective function values under the same control but with different $\tau$ values.

Now, we evaluate the realized reliability level ${\alpha }_{real}$ of the solutions. This can be calculated as a percentage of the resulting temperature realizations that satisfy the temperature constraint under the optimal control sequence $u_{(0,H]}^{*}$ and given samples of $\xi$:

$$\begin{array}{ccc}& & {\alpha }_{real}=\frac{1}{N}\sum _{i=1}^N\mathbb{I}(u^{*}(t,x),{\xi }_i),\hfill \end{array}$$

(43)

$$\begin{array}{ccc}& & \mathbb{I}(u^{*}(t,x),{\xi }_i)=\left\{\begin{array}{cc}1,\hfill & T^{*}(u^{*}(t,x),{\xi }_i,t,x)\le T_{max}\hfill \\ 0,\hfill & T^{*}(u^{*}(t,x),{\xi }_i,t,x)>T_{max}.\hfill \end{array}\right.\hfill \end{array}$$

The result shows that, with $\tau =5\times {10}^{-5}$, we have ${\alpha }_{real}^{OA}={\alpha }_{real}^{IA}=0.973$. It means that, despite the fact that the maximum allowed temperature is very close to the desired temperature, which even in the deterministic case significantly makes it difficult to find admissible values of control, the chance constraint resulted from the optimal control holds.

The computation time for solving NLP problems on the first prediction horizon, for $\tau =5\times {10}^{-5}$, is on average 12.75 h. The decrease of $\tau$ and the corresponding tightening of the feasible set is accompanied by a certain computational penalty, which is presented in

Table 3. Time consumption in hours for the solution of inner and outer problems with decreasing $\tau$.

Problem	$\mathit{\tau }={10}^{-1}$	$\mathit{\tau }={10}^{-2}$	$\mathit{\tau }={10}^{-3}$	$\mathit{\tau }={10}^{-4}$	$\mathit{\tau }=5\times {10}^{-5}$
OA	1.25	3.67	10.31	10.65	10.18
IA	1.14	9.29	15.32	15.04	15.20

. Here, we show only the time taken by the IPOPT optimization and discard the time needed to generate and compile the finite element code. The increase in computation time is possibly also caused by the vanishing gradient problem described above.

From Table 3 and Figure 8, we can see the effect of $\tau$ values in this case. It is seen from Figure 8 that only for $\tau \le {10}^{-3}$ the solutions of inner and outer problems lie in the desired range, i.e., they reach $T_d$ but do not violate $T_{max}$. Table 3 shows that by $\tau \le {10}^{-3}$, the computation needed for solving one problem takes approximately the same time. Thus, the computation time does not heavily depend on the $\tau$ value.

6.2. Results of MPC

Now, we present the solutions of the MPC scheme; i.e., $u_{({\widehat{t}}_1,{\widehat{t}}_1+H^{\prime }]}^{*}(t,x)$ of the inner and outer problems of CCHTMPC, obtained by collecting the first elements of optimal control sequences for all prediction horizons $u_{({\widehat{t}}_z,{\widehat{t}}_z+H]}^{*}(t,x),{\widehat{t}}_z\in ({\widehat{t}}_1,{\widehat{t}}_1+H^{\prime }]$. The number of samples used for the random parameter is $N=300$. The approximation parameter is set to the smallest possible value $\tau =5\times {10}^{-5}$. We decrease it from one prediction horizon to the next one by $\Delta \tau ={10}^{-9}$.

The resulting mean profile of the optimal temperature realizations $E\left[T^{*}(u_{({\widehat{t}}_1,{\widehat{t}}_1+H^{\prime }]}^{*}(t,x),\xi ,t,x)\right]$ over time is shown in Figure A3. It is similar to the temperature profile of the optimal control shown in Figure 7. The difference between the inner and outer solutions is hardly visible; i.e., the mean temperature realizations in the tumor domain converge by the MPC scheme.

To investigate the feasibility of the resulting solutions, the maximum values of the mean temperature distributions on all prediction horizons $E\left[T^{*}(u_{({\widehat{t}}_z,{\widehat{t}}_z+H]}^{*}(t,x),\xi ,t,x)\right],{\widehat{t}}_z\in ({\widehat{t}}_1,{\widehat{t}}_1+H^{\prime }]$ are presented by the gray curves in Figure 10, and the maximum values of the mean temperature distributions of CCHTMPC $E\left[T^{*}(u_{({\widehat{t}}_1,{\widehat{t}}_1+H^{\prime }]}^{*}(t,x),\xi ,t,x)\right]$ (red and blue curves), respectively. It can be clearly seen that the CCHTMPC solution (of the MPC) differs from the CCPDE (of the optimal control). By using MPC, the initial impulse tending to violate the state constraint is not as strong as before, i.e., the temperature does not reach $T_{max}$. At the same time, the solutions after that impulse become much smoother and do not fluctuate anymore.

To evaluate the satisfaction of the chance constraints by the MPC, we plot the temperature realizations of the inner and outer approximate solutions of the CCHTMPC problem for all $N=300$ samples over time; i.e., $T^{*}(u_{({\widehat{t}}_1,{\widehat{t}}_1+H^{\prime }]}^{*}(t,x),{\xi }_i,t,x),i=1,\dots ,N$ with gray curves in Figure 11a for the inner problems and in Figure 11b for the outer problem as well as their mean temperature distributions $E\left[T^{*}(u_{({\widehat{t}}_1,{\widehat{t}}_1+H^{\prime }]}^{*}(t,x),\xi ,t,x)\right]$.

As can be seen, only a small percentage of the trajectories violates the state constraint; i.e., the computed control sequence is feasible for most samples of the random variable. It corresponds to the realized reliability levels ${\alpha }_{real}^{OA,MPC}={\alpha }_{real}^{IA,MPC}=0.987$. Both realized reliability levels are higher than the specified value $\alpha =0.95$ and are also higher than those of the open-loop optimal control ${\alpha }_{real}^{OA}={\alpha }_{real}^{IA}=0.973$. In general, the reliability level should be specified according to the constraints and severity of randomness. The stronger they are, the harder it is to find a feasible solution. In accordance with this, the chosen reliability level must relax the problem enough to have a solution, while providing a sufficient level of performance.

The computation time of CCHTMPC NLP problems using a solution of open-loop optimal control as an initial guess, i.e., on the single prediction horizon for $\tau =5\times {10}^{-5}$ is on average 2.7 h, which is much faster than the computation time of open-loop control. All this demonstrates the significant advantages of MPC due to its feedback feature. However, the computation time still has to be decreased for a real application, as a long computation time results in a long MPC sampling time, which makes the MPC unrealistic.

7. Conclusions and Future Work

The CCMPC of PDE systems poses an enormous computational challenge. In this study, a smoothing inner–outer approximation approach is proposed for the solution of CCMPC of parabolic PDE systems. Through a series of transformations, the original problem is converted to a couple of high-dimensional deterministic state-constrained NLP problems, the solution of which is an approximation to that of the initial problem. Each of these NLPs is solved in a double loop for updates of the initial condition to employ the receding horizon strategy and for decreasing the approximation parameter $\tau$ to guarantee convergence of the approximation to the original problem. The Moreau–Yosida regularization is applied to address the state constraints of the resulting NLPs. The proposed MPC scheme gains the following advantages: (1) feasibility is guaranteed, irrespective of the distribution of input uncertainties (bounded or unbounded uncertainties); (2) through the construction of deterministic tubes, the feasibility is ensured on each prediction horizon, yielding robust constraint satisfaction; and (3) the deterministic tubes are constructed automatically, simplifying derivation of a feasible solution. To verify the properties of the proposed approach, a case study of hyperthermia cancer treatment is performed. Numerical results demonstrate the viability of the proposed method, validate the stated properties, and also indicate some directions for further studies.

The approach turned out to be quite computationally demanding, which points to the need for an accurate temporal discretization and a prediction horizon reduction. The latter can be carried out by analyzing the stability properties and deriving a minimal stabilizing prediction horizon for the MPC. These and other questions, such as reducing the computation time, application-based specification of parameters associated with chance constraints, and discretization improvement, e.g., estimation of errors, and use of spectral techniques will be the focuses of our future research.