![\[ \begin{array}{rl} \min & g(y) + h(u) \\ \mathrm{s.t.} & \left\{\begin{array}{rcl} -\Delta y &=& Bu + f, \mathrm{in\ }\Omega \\ y|_\Gamma &=& y_b\\ u \geq 0 \end{array}\right. \end{array} \]](form_323.png)
Co-Workers:
where
![\[ \begin{array}{rcl} g(y) &=& \displaystyle\int_\Omega (y-y_0)^2 dx \\ h(u) &=& \displaystyle\int_\Omega (u-u_0)^2 dx \end{array} \]](form_324.png)
In the unknowns, 
is control,
is state. This problem has the following optimality condition
![\[ \begin{array}{rcll} -\Delta y &=& Bu + f, & y|_\Gamma = y_b\\ -\Delta p &=& g^\prime(y), & p|_\Gamma = 0 \\ (h^\prime(u) + B^\star p, v-u)&\geq& 0,& \forall v \geq 0 \end{array} \]](form_325.png)
where 
is so-called co-state. It is obvious that the unknowns have very different regularity and often the locations of singularity of control and state are different. Therefore we will adopt different finite element spaces to approximate those unknowns and at the same time, mesh adaptation is implemented in our algorithm to resolve the singularity of the solution. For such approximation, the following posteriori error estimator is obtained.
![\[ \begin{array}{c} \eta_1^2 = \displaystyle\int_{\Omega_0^-} | h^\prime(u) + B^\star p |^2 dx \\ \eta_2^2 = \displaystyle\sum_\tau \displaystyle\int_\tau h_\tau^2 (\mathrm{div} \nabla p + g^\prime (y))^2 dx + \displaystyle\sum_{l\bigcap\Gamma=\emptyset} \displaystyle\int_l h_l [\nabla p \cdot \mathbf{n}]^2 dl \\ \eta_3^2 = \displaystyle\sum_\tau \displaystyle\int_\tau h_\tau^2 (f + Bu + \mathrm{div} \nabla y)^2 dx +\displaystyle\sum_{l\bigcap\Gamma=\emptyset} \displaystyle\int_l h_l [\nabla y \cdot \mathbf{n}]^2 dl \end{array} \]](form_327.png)
The difficulty in the way to approximate the control, the state and the co-state in different finite element spaces is the implementation of the algorithm. We developped the adaptive finite element package AFEPack to solve our problem which provided the ability to handle multi-mesh operation.
Expect for the residual type posteriori error estimator given above, a recovery type posteriori error estimator is also given and it is proved that super-convergence can be obtained the finite element solution of the optimal control problem. The mesh obtained using the recovery type posteriori error estimator is almost the same as the mesh obtained using the residual type posteriori error estimator. Furthmore, the recovery type posteriori error estimator is about to coincide with the numerical error in quantity instead of in profile which is the most desired phenomenon.
In the following, there are some numerical examples.
The exact solution of this example is as
![\[ \begin{array}{rcl} z &=& \left\{\begin{array}{ll} 0.5, & x_1 + x_2 > 1.0 \\ 0.0, & x_1 + x_2 \leq 1.0 \end{array}\right. \\ p &=& \sin \pi x_1 \sin \pi x_2 \\ u_0 &=& 1.0 - \sin {\frac {\pi x_1} 2} - \sin {\frac {\pi x_2} 2} + z \\ u &=& \max(u_0 - p, 0) \\ y_0 &=& 100 \sqrt{(x_1 - 1)^2 + (x_2 - 1)^2} \\ f &=& 4 \pi^4 p - u - 10000/y_0 \\ y &=& 2 \pi^2 p + y_0 \end{array} \]](form_328.png)
and in the following table are the data obtained
![\[ \begin{tabular}{|c|l||c|c|c||c|c|c|} \hline \multicolumn{2}{|c||}{} & \multicolumn{3}{c||}{$y, p$ piecewise linear} & \multicolumn{3}{c|}{$y, p$ piecewise quadratic} \\ \cline{3-8} \multicolumn{2}{|c||}{\raisebox{1.5ex}[0pt]{Item}} & $u$ & $y$ & $p$ & $u$ & $y$ & $p$ \\ \hline & \# nodes & 755 & 972 & 972 & 743 & 245 & 245 \\ \cline{2-8} mesh & \# sides & 1935 & 2725 & 2725 & 1905 & 648 & 648 \\ \cline{2-8} info & \# elements & 1181 & 1754 & 1754 & 1163 & 404 & 404 \\ \cline{2-8} & \# dofs & 3543 & 972 & 972 & 3489 & 921 & 921 \\ \hline \multicolumn{2}{|c||}{$L^2$ error} & 4.88824e-3 & 2.03878e-2 & 1.64006e-3 & 4.88893e-3 & 1.99436e-3 & 1.11350e-4 \\ \hline \end{tabular} \]](form_329.png)
Figures from numerical result is followed
and , color is the value of The next table is the recovery type error estimator and numerical error on sequential adaptive mesh.
![\[ \begin{tabular}{|c|c|c|c|c|c|c|} \hline adaptive step & 1 & 2 & 3 & 4 & 5 \\ \hline \# nodes($u$) & 145 & 204 & 398 & 519 & 721 \\ \hline \# nodes($y$,$p$) & 145 & 537 & 804 & 804 & 804 \\ \hline $||u-u_h||_{L^2}$ & 2.08452e-02 & 1.32012e-02 & 9.07738e-03 & 6.65461e-03 & 5.09988e-03 \\ \hline $|y-y_h|_{H^1}$ & 4.67225e+00 & 2.34528e+00 & 2.07212e+00 & 2.07212e+00 & 2.07212e+00 \\ \hline $|p-p_h|_{H^1}$ & 2.38026e-01 & 1.18983e-01 & 1.05092e-01 & 1.05092e-01 & 1.05092e-01 \\ \hline $||R_hu_h-u_h||_{L^2}$ & 1.79538e-02 & 1.23715e-02 & 8.81438e-03 & 6.11814e-03 & 4.62809e-03 \\ \hline $||G_hy_h-\nabla y_h||_{L^2}$ & 4.89647e+00 & 2.38738e+00 & 2.09266e+00 & 2.09266e+00 & 2.09266e+00 \\ \hline $||G_hp_h-\nabla p_h||_{L^2}$ & 2.45173e-01 & 1.20586e-01 & 1.05759e-01 & 1.05759e-01 & 1.05759e-01 \\ \hline \end{tabular} \]](form_330.png)
where the operator are defined in our paper.
This constraint equation is a weak nonlinear equation as
![\[ \begin{array}{rl} \min & g(y) + h(u) \\ \mathrm{s.t.} & \left\{\begin{array}{rcl} -\Delta y + y^3 &=& Bu + f, \mathrm{in\ }\Omega \\ y|_\Gamma &=& y_b\\ u \geq 0 \end{array}\right. \end{array} \]](form_331.png)
thus the dual equation now looks as
![\[ -\Delta p + 3 y^2 p = g^\prime(y) \]](form_332.png)
The exact solution of our numerical exmaple is as
![\[ \begin{array}{rcl} y_0 &=& \sin 2 \pi x_1 + \sin 2 \pi x_2 \\ y &=& y_0 \\ u_0 &=& \max(4 \pi^2 y_0, 0) \\ u &=& u_0 \\ f &=& 4 \pi^2 y_0 + y_0^3 - u \\ p &=& 0 \\ \end{array} \]](form_333.png)
Data obtained from the numerical result are listed in the following form
![\[ \begin{tabular}{|c|l||c|c|c||c|c|c|} \hline \multicolumn{2}{|c||}{} & \multicolumn{3}{c||}{$y, p$ piecewise linear} & \multicolumn{3}{c|}{$y, p$ piecewise quadratic} \\ \cline{3-8} \multicolumn{2}{|c||}{\raisebox{1.5ex}[0pt]{Item}} & $u$ & $y$ & $p$ & $u$ & $y$ & $p$ \\ \hline & \# nodes & 4473 & 341 & 341 & 4486 & 132 & 132 \\ \cline{2-8} mesh & \# sides & 13006 & 933 & 933 & 13046 & 340 & 340 \\ \cline{2-8} info & \# elements & 8534 & 593 & 593 & 8561 & 209 & 209 \\ \cline{2-8} & \# dofs & 25602 & 341 & 341 & 25683 & 497 & 497 \\ \hline \multicolumn{2}{|c||}{$L^2$ error} & 3.55341e-2 & 1.78840e-2 & 1.98911e-4 & 3.55327e-2 & 3.38543e-3 & 2.02014e-5 \\ \hline \end{tabular} \]](form_334.png)
and , color is the value of The data in the following table are the error estimator and numerical error obtained on sequential adaptive mesh:
![\[ \begin{tabular}{|c|c|c|c|c|c|c|} \hline adaptive step & 1 & 2 & 3 & 4 & 5 \\ \hline \# nodes($u$) & 145 & 384 & 1288 & 2710 & 4290 \\ \hline \# nodes($y$,$p$) & 145 & 174 & 174 & 174 & 174 \\ \hline $||u-u_h||_{L^2}$ & 8.13922e-01 & 2.70932e-01 & 9.71366e-02 & 4.01130e-02 & 2.27846e-02 \\ \hline $|y-y_h|_{H^1}$ & 9.91391e-01 & 9.12057e-01 & 9.12057e-01 & 9.12057e-01 & 9.12057e-01 \\ \hline $|p-p_h|_{H^1}$ & 2.58517e-03 & 2.22193e-03 & 2.22193e-03 & 2.22193e-03 & 2.22193e-03 \\ \hline $||R_hu_h-u_h||_{L^2}$ & 6.95352e-01 & 2.24280e-01 & 8.20353e-02 & 3.09396e-02 & 1.26466e-02 \\ \hline $||G_hy_h-\nabla y_h||_{L^2}$ & 1.03237e+00 & 9.49570e-01 & 9.49570e-01 & 9.49570e-01 & 9.49570e-01 \\ \hline $||G_hp_h-\nabla p_h||_{L^2}$ & 5.20497e-04 & 4.48978e-04 & 4.48978e-04 & 4.48978e-04 & 4.48978e-04 \\ \hline \end{tabular} \]](form_335.png)
The exact solution designed is as
![\[ \begin{array}{rcl} z &=& \left\{\begin{array}{ll} 0.4, & x_1 + x_2 + x_3 > 1.0 \\ 0.0, & x_1 + x_2 + x_3 \leq 1.0 \end{array}\right. \\ p &=& \sin \pi x_1 \sin \pi x_2 \sin \pi x_3 \\ u_0 &=& 1.0 - \sin {\frac {\pi x_1} 2} - \sin {\frac {\pi x_2} 2} - \sin {\frac {\pi x_3} 2} + z \\ u &=& \max(u_0 - p, 0) \\ y_0 &=& 0 \\ f &=& 9 \pi^4 p - u0\\ y &=& 3 \pi^2 p + y_0 \end{array} \]](form_336.png)
The numerical results are in the following table:
![\[ \begin{tabular}{|c|l||c|c||c|c|c|} \hline \multicolumn{2}{|c||}{} & \multicolumn{2}{c||}{On uniform mesh} & \multicolumn{3}{c|}{On adaptive mesh} \\ \cline{3-7} \multicolumn{2}{|c||}{} & Mesh No 1 & Mesh No 2 & Mesh No 1 & Mesh No 2 & Mesh No 3 \\ \hline & \# nodes & 17389 & 132489 & 2315 & 8940 & 31128 \\ \cline{2-7} $u$ & \# edges & 115100 & 900936 & 13080 & 51737 & 182319 \\ \cline{2-7} mesh & \# faces & 192112 & 1523648 & 20083 & 80451 & 285704 \\ \cline{2-7} info & \# elements & 94400 & 755200 & 9317 & 37653 & 134512 \\ \cline{2-7} & \# dofs & 377600 & 3020800 & 37268 & 150612 & 538048 \\ \hline & \# nodes & 17389 & 132489 & 9408 & 12449 & 12724 \\ \cline{2-7} $y,p$& \# edges & 115100 & 900936 & 58351 & 81512 & 83811 \\ \cline{2-7} mesh & \# faces & 192112 & 1523648 & 94150 & 135284 & 139569 \\ \cline{2-7} info & \# elements & 94400 & 755200 & 45206 & 66220 & 68481 \\ \cline{2-7} & \# dofs & 17389 & 132489 & 9408 & 12449 & 12724 \\ \hline & \multicolumn{1}{|c||}{$u$} & 3.15799e-3 & 2.12281e-3 & 2.29716e-3 & 1.87182e-3 & 1.32444e-3 \\ \cline{2-7} $L^2$ error & \multicolumn{1}{|c||}{$y$} & 4.83516e-2 & 1.19974e-2 & 1.28058e-1 & 8.04137e-2 & 7.81510e-2 \\ \cline{2-7} & \multicolumn{1}{|c||}{$p$} & 3.09591e-3 & 7.67661e-4 & 7.77780e-3 & 4.94263e-3 & 4.82434e-3 \\ \hline \end{tabular} \]](form_337.png)
and , color is the value of
![\[ \begin{array}{rl} \min & \displaystyle\int_0^T g(y) + h(u) dt \\ \mathrm{s.t.} & \left\{\begin{array}{rcl} y_t-\Delta y &=& Bu + f, \mathrm{in\ }\Omega \\ y|_\Gamma &=& y_b\\ y|_{t=0} &=& y(0)\\ u \geq 0 \end{array}\right. \end{array} \]](form_338.png)
and its corresponging optimality condition is as
![\[ \begin{array}{rcl} y_t-\Delta y &=& Bu + f \\ -p_t-\Delta p &=& g^\prime(y) \\ (h^\prime(u) + B^\star p, v-u) &\geq & 0, \forall v\geq 0 \end{array} \]](form_339.png)
The character of this problem is in the dual equation, the equation is a parabolic equation in the inverse temporal direction. Now we adopted a layer by layer discrization for the temporal discretization such forward Euler scheme or backward Euler scheme. It is interesting that the forward Euler scheme is almost the same as backward Euler scheme.
Our new point of view to this problem is: The optimal numerical scheme for this problem can be only achieved by a mixed discretization for spatial and temporal dimensions. Our general believings is that if the simulation algorithm is very mature, it will be very useful for optimal control. While according to our new point of view, even if the algorithm for simulation is very stable and efficient, it can be absolutely useless for optimal control at all.
We have got the posteriori error estimator for the Discontinuous Galerkin discritization in temporal direction and for the mixed discretization for spatial and temporal directions. The formation of those estimation are extermely complex.
My research on this topic is partial supported by National Science Funding, China.
![\[ \begin{array}{rl} \min & \displaystyle\frac{1}{2}\displaystyle\int_\Omega (y-y_0)^2 dx + \displaystyle\frac{\epsilon}{2}\displaystyle\int_\Omega u^2 dx \\ \mathrm{s.t.} & \left\{\begin{array}{l} -\Delta y + u y = f, \qquad\mathrm{in\ }\Omega\\ y|_\Gamma = y_b\\ u \geq \delta \end{array}\right. \end{array} \]](form_340.png)
where 
is a negtive constant depended on the coefficient in the Poincare inequality. The first order optimality condition is as
![\[ \begin{array}{rcl} -\Delta y + u y &=& f \\ -\Delta p + u p &=& y-y_0 \\ (\epsilon u - yp, v-u) &\geq & 0,\forall v \geq \delta \end{array} \]](form_342.png)
A posteriori error estimator has been derived using standard finite element discretization. The difficulty arising when 
is a very small number in solving algorithm designment. An efficient enough solving algorithm which is not avaiable yet is our base to give very good numerical result.
1.4.7