Deep multi-input and multi-output operator networks method for optimal control of PDEs

1.   Introduction

The optimal control problem has been successfully applied in various fields, such as heat transfer phenomena ^[1], finance ^[2], image processing ^[3], shape optimization ^[4,5], aerodynamics ^[6,7], crystal growth ^[8], and drug delivery ^[9]. To solve partial differential equation constrained (PDE-constrained) optimal control problems, numerous numerical methods have been developed, including finite element methods, finite differences methods, finite volume methods, spectral methods, and mesh-less methods (see, e.g., ^{[10,11,12,13]}). Despite their effectiveness, the optimal control problem remains challenging to solve, particularly when the problem is nonlinear. Recently, deep learning has emerged as a popular method for solving partial differential equations, especially nonlinear ones.

Deep learning methods for solving PDEs have received significant attention, including physics-informed neural networks (PINN) ^[14,15], deep Galerkin method ^[16], deep Ritz method ^[17], deep Nitsche method ^[18], and deep operator networks (DeepONets) ^[19,20]. PINNs can be used to solve specific PDEs with boundary conditions and loading terms, but they require expensive optimization computation cost during inference. Therefore, the PDEs with operating conditions and real-time inference cannot be solved by PINNs. The neural operator full-fills well. Various versions of operator networks have been published, e.g., graph neural operator networks ^[21], Fourier neural operator (FNO) ^[22], physics-informed neural operators (PINO) ^[23], and deep multiple input operator network (DeepMIONet) ^[24].

Neural networks have several advantages over traditional numerical solvers, including being mesh-free and easier to deal with complex geometric regions. Simulating control problems that involve these complex geometric regions using traditional numerical methods often requires high-quality grids and extensive preprocessing before simulation. By contrast, researchers are interested in using neural networks to replace traditional numerical methods. To solve optimal control problems with PINNs, methodologies and guidelines have been proposed in previous works ^[25,26]. Barry-Straume et al. used a two-stage framework to solve PDE-constrained optimization problems ^[27]. Wang et al. used physics-informed deep operator networks (DeepONets) to learn the solution operator of parametric PDEs, building a surrogate for solving PDE-constrained optimization problems ^[28]. In summary, deep learning approaches, such as PINNs and DeepONets, have shown promise for solving PDE-constrained optimal control problems, providing an efficient solution without requiring extensive preprocessing or expensive optimization during inference. Future work could explore further improvements to deep learning methods and investigate their application in new fields.

In this paper, to solve PDE-constrained optimal control problems with available data, we introduce a MIMOONet. In the context of the PDE optimal problem, the governing PDE is fully known, and the objective is to determine a control variable that minimizes the cost function. Initially, the PDE-constrained optimal control problem is reformulated into a PDE system using the adjoint method. Subsequently, the PDE system is tackled by using MIMOONets. Additionally, we examine a physical system described by PDEs and propose physics-informed MIMOONets for addressing the PDE-constrained optimal control problem. Overall, this method (MIMOONet) has the following advantages:

● MIMOONet can solve PDE-constraint problems governed by different types of PDE.

● MIMOONet can easily approximate nonlinear optimal control problems.

● Compared with traditional numerical methods, its prediction speed is faster.

The remainder of this paper is organized as follows. In Section 2, we introduce the adjoint state method of solving PDE-constrained optimization problem and optimality system, and provide the framework of MIMOONets, physics-informed MIMOONets, and the detailed method of our main technical contribution. In Section 3, we give the deep learning framework of elliptic and parabolic constrained optimal control problem, and present the numerical results to assess the performance of the proposed MIMOONets and physics-informed MIMOONets. Finally, Section 4 summarizes the results, potential pitfalls, and shortcomings, and details the groundwork for future directions.

2.   Methods

Let $U$, $S$, $V$ be Banach spaces. We consider the following PDE-constrained optimization problem:

where $J: S\times U \rightarrow \mathbb{R}$ is a cost function, $\mathcal{F}: S\times U \rightarrow V$ is a system of PDEs subject to initial and boundary conditions, $u$ and $v$ are the state variable and control variable, respectively. Assume that the problem (2.1) has unique solutions $u$ and $v$. In the subsequent sections, the MIMOONets method is presented under this assumption.

2.1.   Optimality system

In this subsection, we transform the PDE-constrained optimization problem into an optimality system.

The PDE-constrained optimization problem, which involves optimizing an objective function subject to a set of partial differential equations (PDEs), can be transformed into an optimality system. The resulting optimality system consists of two sets of equations: the state equations and the adjoint equations. These equations are coupled and must be solved simultaneously to obtain the optimal solution to the original PDE-constrained optimization problem.

Consider the following problem ^[29],

where $x$ and $t$ denote space and time variables, respectively, the domain $\Omega\subseteq\mathbb{R}^d$, $\partial\Omega$ is the boundary of the domain $\Omega$, and $\mathcal{B}$ and $\mathcal{I}$ are boundary conditions and initial condition, respectively. We construct the Lagrangian function for the problems (2.2) and (2.3) as follows:

Here, $p_{1}, p_{2}$, and $p_{3}$ are Lagrange multiplier functions defined on $\Omega\; \times\; [0, T]$, $\partial\Omega\; \times\; [0, T]$ and $\partial\Omega\; \times\; {0}$, respectively. According to the Lagrange principle, we seek the pair $(u^{*}, v^{*})$ and the Lagrange multipliers or adjoint field $\mathbf{p} = (p_{1}, p_{2}$, $p_{3})$ to satisfy the optimality conditions. Therefore, the problems (2.2) and (2.3) are equivalent to the following unconstrained problem

Then, the directional derivative of $\mathcal{L}$ with respect to $u$ disappears at the optimal point, that is

For the control variable $v$ and Lagrange multipliers $\mathbf{p}$, we have

and

Therefore, the problems (2.2) and (2.3) can be written in the following way:

Once the optimality system has been derived, it can be solved by using neural networks. Next, we introduce the MIMOONet and the physics-informed MIMOONet methods for solving the optimality system.

2.2.   Multiple-input operators networks

DeepONets is a learning framework proposed by Lu et al. ^[19], which enables the learning of abstract nonlinear operators in infinite-dimensional function spaces. It is inspired by the universal approximation theorem of operators ^[30]. The DeepONets network comprises two main components: the trunk network and the branch network.

The trunk network provides the basis functions for the output function by encoding information related to the space-time coordinates. It takes as input the space-time coordinates and any other relevant physical parameters and produces a set of basis functions for the function space. These basis functions serve as a representation of the output function and are used in the computation of the final output.

The branch network encodes the input function to provide the coefficients at fixed sensor points. Given an input function, which may be a solution to a PDE or a data-driven function, the branch network maps it to a set of coefficients that correspond to specific sensor points in the domain. These coefficients are then combined with the basis functions from the trunk network to produce the final output function.

Theorem 2.1. Suppose that $X$ is a Banach space, $K_{1} \subset X$, $K_{2}\subset R^{d}$ are two compact sets in $X$ and $R^{d}$, respectively. Let $V$ be a compact set in $C(K_{1})$, $G:V\rightarrow C(K_{2})$ a nonlinear continuous operator, $\sigma$ a continuous non-polynomial function. Then, for $\forall\varepsilon > 0$, there exist $n$, $p$, $m\in \mathbb{N}$, constants $c^{k}_{i}$, $a^{k}_{ij}$, $\theta_{i}^{k}$, $\zeta_{k}\in R$, $w_{k}\subset R^{d}$, $x_{j}\in K_{1}$, $i = 1, 2, \cdots, n$, $k = 1, \cdots, p$, $j = 1, \cdots, m$, such that

for any $u\in V$ and $y\in K_{2}$.

DeepONets are designed to learn abstract nonlinear operators in infinite-dimensional function spaces from a single input function defined on a Banach space, while real-world problems often involve multiple input functions. To address this issue, the DeepMIONet was proposed in ^[24], which is defined through the tensor product of Banach spaces.

In DeepMIONets, only the case of two input functions $u$ and $v$ are considered, which are represented by their respective branch networks. The trunk network provides basis functions that span the function space, and the outputs of the branch networks are combined with the basis functions by using tensor products to produce the final output function. Specifically, the output of DeepMIONet $G_{\theta}(u, v)$ can be written as:

where $b^{u}_{i}$ and $b^{v}_{i}$ denote the $i$-th output of the branch networks corresponding to the input functions represented by $u$ and $v$, respectively. And $tr_{i}$ is the $i$-th output of the trunk network.

The architecture of DeepMIONets includes two separate branch networks and a shared trunk network (as shown in Figure 1), which provide basis functions spanning the function space. The branch networks encode the input functions and provide coefficients at fixed sensor points, while the trunk network provides basis functions spanning the function space. The outputs of the branch networks and the trunk network are combined using tensor products to produce the final output function. Specifically, the outputs of the branch networks are the coefficients of the input functions, while the trunk network provides a set of basis functions that are combined with the outputs of the branch networks using tensor products to obtain the final output function.

2.3.   Multi-input and multi-output operators networks

Here, we focus on using neural operator networks to solve systems of PDEs.

Although DeepONets and DeepMIONet can be used to solve a single PDE, they cannot do it for a system of PDEs, in which at least two solution operators must be generated and two output operators are needed for a network. To solve this problem, we propose MIMOONets, which also can be used to solve systems of PDEs. The MIMOONets framework is composed of a trunk network and multiple branch networks. The trunk network provides the basis functions of the solution operators, while the branch networks provide additional groups of coefficients of the solution operators at fixed sensor points.

We consider the following problem that is a system of PDEs in domain $D\subseteq \mathbb{R}^{d}$ ($d$ is the dimension of space),

where $u_{i}$ and $f_{i}$ are functions, $\mathcal{L}_{i}$ is the differential operator, $\varphi_{i}$ is the boundary condition of $u_{i}$. Let $G^{i}$ be the solution of the operator with input functions $f_{i}$ and $\varphi_{i}$ s.t.

and

This means that $G^{i}(f_{1}, f_{2}, \cdot\cdot\cdot, f_{n}, \varphi_{1}, \varphi_{2}, \cdot\cdot\cdot, \varphi_{n})(y)$ is the corresponding output function. According to the results in ^[24] and ^[30], the solution $u_{i}$ can be expressed as

Then, the solution of problem (2.12) can be learned by using MIMOONets, which is defined through the tensor product of Banach spaces.

Here, we only give the framework of two-input and two-output operators networks (see Figure 2). The solution operator networks can be described by

where the definitions of $b1^{f}_{k}$, $b2^{f}_{k}$, $b1^{g}_{k}$, $b2^{g}_{k}$ and $tr_{k}$ are similar to (2.11). To reduce the generalized error, we may add a bias $b^{1}_{0}, b^{2}_{0}\in\mathbb{R}$ in the last stage:

Let $G^{i}:C(D)\longrightarrow L^{2}(D)$ be a Borel measurable mapping with $G^{i}\in L^{2}(\mu)$. Then, for any $\varepsilon > 0$, there exists an operator network $G^{i}_{\theta}:C(D)\longrightarrow L^{2}(D)$, such that

here, $\mu$ is a probability measure on $C(D)$.

When we have some dataset of the pair solution $\{u(x_{k}), v(x_{k})\}_{1}^{N}$, the solution can be learned by MIMOONets. The corresponding loss function can be formulated as follows

Now we can learn parameter $\theta$ by minimizing the loss function (2.15) and using the stochastic gradient descent method.

2.4.   Physics-informed MIMOONets

A large amount of paired datasets is required to solve PDE systems by using MIMOONets. However, data acquisition is expensive in many engineering applications and physical systems. Under the condition of sparse data, it becomes important to introduce physics-informed neural networks to train MIMOONets by integrating known differential equations with label data in the loss function. We use automatic differentiation of the outputs of MIMOONets with respect to their input coordinates and adopt an appropriate regularization mechanism to ensure that the target output functions satisfy the PDE constraints.

For simplicity, we consider the following problem (without causing confusion, we still use the preceding symbols):

The solutions $u$ and $v$ can be expressed as

For the problem (2.16), the loss function of a physics-informed MIMOONets is defined as follows

where

and

Here, $N$ is the number of initial data points. $N_{f}$ and $N_{g}$ are the number of sample from the computational domain $\Omega$ for the PDEs. $N_{\varphi}$ and $N_{\psi}$ are the number of the boundary points for $u$ and $v$.

The loss function (2.18) is minimized by learning the parameters $\theta$ of the deep neural network. Sometimes, in order to improve the accuracy of the numerical solution or increase the convergence rate, we can apply penalty parameters.

2.5.   PDE-constrained optimization with physics-informed MIMOONets

For a given PDE-constrained optimal control problem such as (2.2), we use physics-informed MIMOONets to solve the optimization problems. The corresponding steps are given as follows. First, we turn the PDE-constrained optimization problem (2.2) to an optimality system (2.9), which consists of adjoint equation, state equation, and optimality condition; second, we solve the optimality system using physics-informed MIMOONets. The detailed computing framework can be seen in Algorithm 1.

Remark 2.1. For MIMOONets (no physics-informed), the algorithm is only need to change "Input: $\cdots$" to "Input: $N$ (the number of initial data of $\{(x_i, u(x_i), v(x_i))\})$".

3.   Main results

In the following demonstrations, to showcase the effectiveness of MIMOONs, some numerical examples of elliptic, semi-linear elliptic, parabolic, and second-order hyperbolic optimal control problems are provided. Data-driven MIMOONs or physics-informed MIMOONs are employed with uniform distribution random sampling on the solution domain during the training process. In the subsequent examples, the relative error of the numerical solution $\boldsymbol{u}$ is calculated by: $\dfrac{||\boldsymbol{u}- \boldsymbol{u}^{\ast}||}{||\boldsymbol{u}^{\ast}||}$, where the reference solution $\boldsymbol{u}^{\ast}$ is an analytical solution, or finite element approximated solution with $100\times100$ spatial grid.

3.1.   Linear elliptic optimal control problem

We start with an example involving finding an optimal heat source under homogeneous Dirichlet boundary conditions. The model can be represented as follows,

where $\Omega$ is a bounded domain, $u:\Omega\rightarrow \mathbb{R}$ is the unknown temperature satisfying (3.2), $u_{d}:\Omega\rightarrow \mathbb{R}$ is the given desired temperature, $v$ is the unknown control function, $f$ is the source term in $\Omega$, $\alpha\geq 0$ is a regularization parameter. Here, we set $\Omega = (0, 1)\times (0, 1)$, $u_{d}(x, y) = (1-10\pi)\sin(\pi x)\sin(\pi y)$, $\lambda = 1$, $f(x, y) = (5+2\pi^2)\sin(\pi x)\sin(\pi y)$. When, $u(x, y) = \sin(\pi x)\sin(\pi y)$, $v(x, y) = -5\sin(\pi x)\sin(\pi y)$, the $J(u, v)$ gets the global minimum. We know that the optimal control problems (3.1) and (3.2) can be transformed into the optimality system as follows

We use the MIMOONets to solve the PDEs system (3.3). The solution operator $G^{1}$ and $G^{2}$ of $f$ and $u_{d}$ can be represented as follows,

where the branch network 1 and the branch network 2 are two separate 5-layer fully connected neural networks (FNN). Every hidden layer and output layer has 300 neurons. The input layer per network has 100 neurons. The trunk network is 5-layer FNN with 300 neurons per hidden layer and 150 neurons output layer. Relu or Tanh is used as the activation function. The loss function of the deep MIMOONets is denoted by

We use the same network structure as MIMOONets to consider the physics-informed MIMOONets. The activation function is Tanh. The corresponding loss function can be expressed as follows

where

and

where $N$ is the number of initial data points. $N_{f}$ and $N_{u_{d}}$ denotes the number of integration points of $f$ and $u_{d}$ in the computational domain $\Omega$, respectively. $N_{BC}$ is the number of the boundary points on the $\partial\Omega$.

To evaluate the loss, we randomly sample $N = 10,000$ training points $(x_{k}, y_{k})\in \Omega$, $N_{f} = N_{u_{d}} = 10,000$ residual training points $(x_{k}, y_{k})\in \Omega$. We select $N_{BC} = 400$ equidistant boundary training points $(x_{i}, y_{i})\in \partial\Omega$.

We utilize the Adam optimizer from PyTorch to train the networks, and take the corresponding learning rate 0.001. At the same time, we can choose different parameters to improve the accuracy of deep physics-informed MIMOONets method for the following format

We take $\lambda_{0} = \lambda_{1} = \lambda_{2} = 1$, $\lambda_{3} = \lambda_{4} = 100$. The experimental results are shown in Table 2 and Figure 5.

Figures 3–5 show that both the deep MIMOONets method and the deep physics-informed MIMOONets method are effective for optimal control problems with elliptic constraints; also, the activation function Relu converges faster than Tanh with the MIMOONets method. However, the precision of the deep MIMOONets method can be achieved with only a small amount of data by means of the deep physics-informed MIMOONets method.

3.2.   Semi-linear elliptic optimal control problem

A semi-linear elliptic constrained optimal control problem is considered here. The model is described as follows

The problems (3.9) and (3.10) lead to the following optimality system

We use a multi-output operators network (MOONet) to solve the PDEs system (3.11). The corresponding solution operators $G_{1}$ and $G_{2}$ are shown below

where the branch network is a 5-layer FNN with 300 neurons per hidden layer and output layer. The input layer contains 100 neurons. The trunk network is composed of five hidden layers with 300 neurons in each layer and 150 neurons in the output layer. The loss function for known data is similar to (3.5). Takeing into account the physics-informed MOONet and $\alpha = 1$, the corresponding loss function is expressed as follows

where

and

When we use physics-informed MOONet, we choose $N = N_{1} = N_{2} = 10,000$, $N_{BC} = 400$, $\lambda_{0} = \lambda_{1} = \lambda_{2} = 1$, $\lambda_{3} = \lambda_{4} = 300$. Not using physics-informed MOONet, we take $N = 10,000$, $\lambda_{0} = 1$, the others are 0. The sample data are collected by finite difference and sequential quadratic programming. The experimental results are shown in Tables 7 and 8.

We find that the Relu of the activation function converges faster than Tanh using the MOONETS method, which is similar to the linear elliptic optimal control problem. The deep physical information MOONets method requires very little data to achieve the same precision of the deep MOONets method.

3.3.   Parabolic optimal control problem

We consider the following parabolic optimal control problem,

where $\Omega = (0, 1)$, $D = \Omega\times(0, T]$, $u:D\rightarrow \mathbb{R}$ is the unknown term satisfying (3.17), $u_{d}:D\rightarrow \mathbb{R}$ is the given desired temperature distribution, $v$ is the unknown control function, $f$ is the source term in $D$, and $\alpha\geq 0$ is a regularization parameter. Here, we set $u_{d}(x, t) = \sin(\pi x)(t+1)+\frac{1}{2}e^{t}\sin(\pi x)-\frac{1}{2}(e^{t}-e)\sin(\pi x)$, $\alpha = 1$, $f(x, t) = \sin(\pi x)+\pi^{2}(t+1)(\sin(\pi x)+\frac{1}{2}(e^{t}-e)\sin(\pi x)$, $\varphi(x) = sin(\pi x)$. When $v(x, t) = -\frac{1}{2}(e^{t}-e)\sin(\pi x)$, $u(x, t) = (t+1)\sin(\pi x)$, the $J(u, v)$ gets the global minimum.

The optimal control problems (3.16) and (3.17) can be transformed into the optimality system,

For the PDEs system (3.18), we use MIMOONets to solve it. The operator $G^{1}$ and $G^{2}$ can be learned from the source term $f$, $u_{d}$ and $\varphi$. Their representations are as follows

Here, we select 100 sensors points for input functions $f$, $u_{d}$ and $\varphi$. Three branch networks are three separate 5-layer FNN with 300 neurons per hidden layer and output layer. The trunk network is a 5-layer FNN with 300 neurons per hidden layer and 150 neurons for output layer. Relu or Tanh is used as the activation function.

We use the same network structure of the MIMOONets as the network structure physics-informed MIMOONets where the activation function is Tanh. The corresponding loss function for the deep MIMOONets is expressed by

The deep physics-informed MIMOONets loss function is the following form

where

and

To calculate the value of the loss function, we take the random sample of initial data points $N = 10,000$, the numbers of the residual training points $N_{f} = N_{u_{d}} = 10,000$, the number $N_{IC}$ and $N_{TC}$ of the initial and termination condition points $x\in\Omega$ are 100, and the number $N_{BC}$ of the boundary condition is 100. Then, we use the Adam optimizer to train the deep MIMOONets and physics-informed MIMOONets ($\lambda_{i} = 1, i = 0, 1, 2, 3;\lambda_{j} = 100, j = 4, 5, 6$) by minimizing the loss of Eqs (3.20) and (3.21). The learning rate is 0.002. The experimental results are shown in Tables 5 and 6, and Figures 6–8.

We found that both the deep MIMOONets method and the deep physics-informed MIMOONets method are effective for parabolic optimal control problem, and the activation function Relu converges faster than Tanh with the MIMOONets method. We also found that the physics-informed MIMOONets not only attains comparable accuracy to the original MIMOONets but also satisfies the underlying PDEs constraint.

3.4.   Hyperbolic optimal control problem

Consider the following hyperbolic optimal control problem ^[13]:

where $\Omega = [0, 1]^2$, $D = \Omega\times(0, T]$. In the experiment, we take $u_{d}(x, t) = \sin(\pi x_{1})\sin(\pi x_{2})(e^t+2+2\pi^2(t-T)^2)$, $T = 1, \alpha = 1, a = 0.2, b = 0.8$, $f(x, t) = (1+2\pi^2)e^t \sin(\pi x_{1})\sin(\pi x_{2})-\max\{a, \min\{b, (t-T)^2\sin(\pi x_{1})\sin(\pi x_{2})\}\}$, $\varphi(x) = \sin(\pi x_{1})\sin(\pi x_{2})$, $\psi(x) = \sin(\pi x_{1})\sin(\pi x_{2})$. The exact solution $u(x, t) = e^t\sin(\pi x_{1})\sin(\pi x_{2})$, $v(x, t) = \max\{a, \min\{b, (t-T)^2\sin(\pi x_{1})\sin(\pi x_{2})\}\}$.

Based on (3.23) and (3.24), the following optimality system can be obtained,

The solution operator of system (3.25) can be represented by $G^{1}$ and $G^{2}$ as follows:

In the experiment, we use 100 sensor points as input functions $f$, $u_{d}$, $\varphi$, and $\psi$. Each of the four branch networks consists of an independent 5-layer FNN, with 300 neurons in each hidden layer and output layer. The trunk network is a 5-layer FNN with 300 neurons per hidden layer and 150 neurons for the output layer. The activation function used is either Relu or Tanh.

We employ the same network structure as the physics-informed MIMOONets, using Tanh as the activation function. The corresponding loss function for the deep MIMOONets is expressed as follows:

The deep physics-informed MIMOONets loss function is the following form

where

and

In the experiment, uniformly random samples are taken with $N = 10,000$ initial data points, and the number of residual training points is $N_{f} = N_{u_{d}} = 125,000$. The sample size for initial condition points $N_{IC}$ and termination condition points $N_{TC}$ in $\Omega$ is 10,000, and the number of points on the spatiotemporal boundary $\partial\Omega\times [0, T]$ is $N_{BC} = 400\times100$. For the sampling of sensor points for input functions, besides the source function $f(x_{1}, x_{2}, t)$, 100 sensor points are uniformly randomly sampled. However, for the source function $f(x_{1}, x_{2}, t)$ in the spatiotemporal space $D$, Latin hypercube sampling is used to sample 100 sensor points. Subsequently, we use the Adam optimizer to train the deep MIMOONets and physics-informed MIMOONets ($\lambda_{0} = 0, \lambda_{i} = 1, i = 1, 2, 3;\lambda_{j} = 100, j = 4, 5, 6$) by minimizing the loss of Eqs (3.27) and (3.28). The learning rate is set to 0.0001. The experimental results are shown in Tables 7 and 8, and Figures 9–11.

For hyperbolic optimal control problems, both the deep MIMOONets method and the deep physics-informed MIMOONets method are effective. In the experiments, the MIMOONets method adopts a data-driven approach, and the Relu activation function converges faster than Tanh. As for the physics-informed MIMOONets method, we adopt a non-data-driven mode and conduct experiments by considering the PDE as a constraint. The experimental results are excellent.

4.   Conclusions

In this work, a novel deep learning framework is presented, which enables the construction of fast surrogates for solving PDE-constrained optimization problems using MIMOONets and physics-informed MIMOONets. The MIMOONets frameworks (MIMOONets and physics-informed MIMOONets) offer flexibility and faster implementation compared to other traditional methods. Compared with MIMOONets, the physics-informed MIMOONets require little paired input-output data, and are more efficient and cost-effective.

Although our methods (MIMOONets and physics-informed MIMOONets) were initially designed for solving PDE-constrained optimization problems, they can also be extended to multi-equation-coupled problems, including hyperbolic and parabolic optimal control problems ^[31], Cahn-Hilliard-Navier-Stokes equation ^[32], and other scenarios. Moreover, they can be effectively employed for real-time prediction in optimal systems governed by PDEs.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work is supported by the National Natural Science Foundation of China (Granted No. 11961008) and Funding from the Scientific Research Fund Project of Guizhou Education University (Granted No. 2024ZD007).

Conflict of interest

The authors declare there is no conflict of interest.

Appendix

A.   Supplementary visualizations of elliptic optimal control problem

We present some numerical images for solving elliptic optimal control problems (3.1) and (3.2) using deep MIMOONets and physics-informed MIMOONets framework. The data-driven results are shown in Figures A.1–A.6. When there is no data, we can also use deep physics-informed MIMOONets to solve the problem. The experimental results are shown in Figure A.7.

B.   Supplementary visualizations of parabolic optimal control problem

We present some numerical images for solving parabolic optimal control problems using deep MIMOONets and physically informed MIMOONets frameworks. The data-driven results are shown in Figures B.1–B.6. When there is no initial data, we can also use MIMOONets based on depth physical information to solve this problem, and the experimental results are shown in Figure B.7.