---
abstract: 'We consider a one dimensional elliptic distributed optimal control problem with pointwise constraints on the derivative of the state. By exploiting the variational inequality satisfied by the derivative of the optimal state, we obtain higher regularity for the optimal state under appropriate assumptions on the data. We also solve the optimal control problem as a fourth order variational inequality by a $C^1$ finite element method, and present the error analysis together with numerical results.'
address:
- 'Susanne C. Brenner, Department of Mathematics and Center for Computation and Technology, Louisiana State University, Baton Rouge, LA 70803, USA'
- 'Li-yeng Sung, Department of Mathematics and Center for Computation and Technology, Louisiana State University, Baton Rouge, LA 70803, USA'
- 'Winnifried Wollner, Department of Mathematics, Technische Universität Darmstadt, 64293 Darmstadt, Germany'
author:
- 'S.C. Brenner'
- 'L.-Y. Sung'
- 'W. Wollner'
title: A One Dimensional Elliptic Distributed Optimal Control Problem with Pointwise Derivative Constraints
---

Ø

[^1]

Introduction {#sec:Introduction}
============

Let $I$ be the interval $(-1,1)$ and the function $J:L_2(I)\times L_2(I)\longrightarrow \R$ be defined by $$\label{eq:JDef}
  J(y,u)=\frac12 \big(\|y-y_d\|_\LT^2+\beta\|u\|_\LT^2\big),$$ where $y_d\in L_2(I)$ and $\beta$ is a positive constant.

The optimal control problem is to $$\label{eq:OCP}
  \text{find}\quad(\bar y,\bar u)=\argmin_{(y,u)\in \bbK}J(y,u),$$ where $(y,u)\in H^1_0(I)\times L_2(I)$ belongs to $\bbK$ if and only if $$\begin{aligned}
{3}
  \int_I y'z'dx&=\int_I (u+f)z\,dx&\qquad&\forall\,z\in H^1_0(I),
  \label{eq:PDEConstraint}\\
  y'&\leq \psi &\qquad&\text{a.e. on $I$}.\label{eq:DerivativeConstraint}\end{aligned}$$ We assume that $$\label{eq:DataRegularity}
  f\in H^1(I), \; \psi\in H^2(I)$$ and $$\label{eq:psiConstraint}
  \int_I \psi\,dx>0.$$

\[rem:Literature\] The optimal control problem defined by – is a one dimensional analog of the optimal control problems considered in [@CB:1988:Gradient; @CF:1993:Gradient; @DGH:2009:Gradient; @OW:2011:Gradient; @Wollner:2012:Gradient]. It was solved by a $C^1$ finite element method in [@BSW:2020:OneD] under the assumptions that $$\label{eq:OldData}
  f\in H^{\frac12-\epsilon}(I)\quad\text{and}\quad
  \psi\in H^{\frac32-\epsilon}(I).$$

Since the constraint implies $y\in H^2(I)$ by elliptic regularity, we can reformulate the optimization problem – as follows: $$\label{eq:ROCP}
  \text{Find}\quad \bar y=\argmin_{y\in K}\frac12\big(\|y-y_d\|_\LT^2
  +\beta\|y''+f\|_\LT^2\big),$$ where $$\label{eq:KDef}
  K=\{y\in H^2(I)\cap H^1_0(I):\,y'\leq \psi \;\text{on}\;I\}.$$

According to the standard theory [@ET:1999:Convex; @KS:1980:VarInequalities], the minimization problem defined by - has a unique solution characterized by the fourth order variational inequality $$\beta\int_I (\bar y''+f)(y''-\bar y'')dx+ \int_I (\bar y-y_d)
   (y-\bar y)dx\geq 0
    \qquad\forall\,y\in K,$$ which can also be written as $$\label{eq:VI}
  a(\bar y,y-\bar y)\geq \int_I y_d(y-\bar y)dx-\beta\int_I f(y''-\bar y'')dx
  \qquad\forall\,y\in K,$$ where $$\label{eq:aDef}
  a(y,z)=\beta\int_I y'z'dx+\int_I yz\,dx.$$

\[rem:VI\] The reformulation of state constraint optimal control problems as fourth order variational inequalities was discussed in [@PS:1996:OC], and a nonconforming finite element based on this idea was introduced in [@LGY:2009:Control]. Other finite element methods can be found in [@GY:2011:State; @BSZ:2013:OptimalControl; @BDS:2014:PUMOC; @BOPPSS:2016:OC3D; @BGPS:2018:Morley; @BGS:2018:Nonconvex; @BSZ:2019:Neumann].

\[rem:Nontrivial\]Note that implies $$\int_I\psi\,dx\geq \int_I y'\,dx=0 \qquad \forall\,y\in K$$ and hence $\int_I \psi\,dx\geq0$ is a necessary condition for $K$ to be nonempty. It is also a sufficient condition because the function $y$ defined by $$y(x)=\int_{-1}^x (\psi(t)-\bar\psi)dt$$ belongs to $K$, where $\bar\psi$ is the mean of $\psi$ over $I$. Furthermore, $$0=\int_I \psi\,dx=\int_I(\psi-y')dx$$ together with implies $\psi=y'$ identically on $I$ and hence $K=\{\psi\}$ is a singleton. Therefore we impose the condition to ensure that the optimization problem defined by – is nontrivial.

Our goal is to show that $\bar y\in H^3(I)$ under the assumptions in and consequently / can be solved by a $C^1$ finite element method with $O(h)$ convergence in the energy norm. Note that previously $\bar y\in H^{\frac52-\epsilon}$ was the best regularity result in the literature for Dirichlet elliptic distributed optimal control problems on smooth/convex domains with pointwise constraints on the gradient of the state.

The rest of the paper is organized as follows. The $H^3$ regularity of $\bar y$ is obtained in Section \[sec:DerivativeVI\] through a variational inequality for $\bar y'$ that can be interpreted as a Neumann obstacle problem for the Laplace operator. The $C^1$ finite element method for / is analyzed in Section \[sec:Discrete\], followed by numerical results in Section \[sec:Numerics\]. We end with some remarks on the extension to higher dimensions in Section \[sec:Conclusions\].

A Variational Inequality for $\bar y'$ {#sec:DerivativeVI}
======================================

Observe that the set $\{y':\,y\in K\}$ is the subset $\cK$ of $H^1(I)$ given by $$\begin{aligned}
\label{eq:cKDef}
  \cK&=\{v\in H^1(I):\,\int_I v\,dx=0\quad\text{and}\quad v\leq\psi
  \;\text{on}\;I\},
      %=\{v\in H^1(I)/\R:\,v\leq\psi \;\text{on}\;I\},\end{aligned}$$ and the variational inequality is equivalent to $$\begin{aligned}
\label{eq:NeumannVI}
  &\int_I (\Phi-f')(q-p)dx+\int_I p'(q'-p')dx\\
    &\hspace{40pt}+[f(1)(q(1)-p(1))-f(-1)(q(-1)-p(-1))]\geq0 \qquad
    \forall\,q\in \cK,\notag\end{aligned}$$ where $p=\bar y'$, $q=y'$, and $\Phi\in H^1(I)$ is determined by $$\begin{aligned}
 \beta\Phi'&=y_d-\bar y \label{eq:PhiDef1}
\intertext{and}
 \int_I\Phi\,dx&=0.\label{eq:PhiDef2}\end{aligned}$$ Moreover is the variational inequality that characterizes the solution of the following minimization problem: $$\label{eq:NeumannObstacle}
  \text{Find}\quad p=\argmin_{q\in \cK}\Big[\frac12 \int_I (q')^2dx
  +\int_I(\Phi-f')q\,dx
   +f(1) q(1)- f(-1) q(-1)\Big].$$

A Neumann Obstacle Problem {#subsec:Neumann}
--------------------------

The minimization problem , which is a Neumann obstacle problem, can be written more conveniently as $$\label{eq:ObstacleProblem}
  p=\argmin_{q\in \cK}\Big[\frac12 b(q,q)+(\phi,q)+ \tau q(1)
  -\sigma q(-1)\Big],$$ where $\sigma=f(-1)$, $\tau=f(1)$, $$\label{eq:Shorthands}
 b(q,r)=\int_I  q'r'\,dx,\quad
 (\phi,q)=\int_I \phi\, q\,dx \quad\text{and}\quad
 \phi=\Phi-f'.$$ Note that we have a compatibility condition $$\label{eq:Compatibility}
 \int_I \phi\,dx+\tau-\sigma=0$$ that follows from , and the Fundamental Theorem of Calculus for absolutely continuous functions.

Since $b(\cdot,\cdot)$ is coercive on $H^1(I)/\R$, the obstacle problem defined by and has a unique solution $p$ characterized by the variational inequality $$\label{eq:VI1}
 b(p,q-p)+(\phi,q-p)+\tau(q(1)-p(1))-\sigma(q(-1)-p(-1))\geq 0
 \qquad\forall\,q\in \cK.$$

\[thm:pRegularity\] The solution $p=\bar y'\in \cK$ of / belongs to $H^2(I)$.

We begin by observing that $$\label{eq:VI2}
 b(p,q-p)+(\phi,q-p)+\tau(q(1)-p(1))-\sigma(q(-1)-p(-1))\geq 0
  \qquad\forall\,q\in\tilde K,$$ where $$\label{eq:tKDef}
  \tilde K=\{q\in H^1(I):\,q\leq\psi\;\;\text{in}\;I\;\;
  \text{and}\;\int_I q\,dx\geq0\}.$$ Indeed, $q\in \tilde K$ implies $q-\bar q\in K$, where $\bar q$ is the mean of $q$ over $I$, and hence, in view of , the definition of $b(\cdot,\cdot)$ in and , $$\begin{aligned}
    &b(p,q-p)+(\phi,q-p)+\tau(q(1)-p(1))-\sigma((q(-1)-p(-1))\\
      &\hspace{10pt}=b(p,q-\bar q-p)+(\phi,q-\bar q-p)+
      \tau(q(1)-\bar q-p(1))-\sigma(q(-1)-\bar q-p(-1))\\
      &\hspace{10pt}\geq 0\end{aligned}$$ for all $q\in\tilde{K}$.

Let $\fK\subset H^1(I)$ be defined by $$\label{eq:bKDef}
 \fK=\{q\in H^1(I):\;q\leq \;\psi\;\;\text{in}\;I\},$$ and $G:H^1(I)\longrightarrow [0,\infty)$ be defined by $$\label{eq:gDef}
   G(q)=\int_I q\,dx.$$ Then the function $\psi$ belongs to $\fK$ and $$\label{eq:Slater}
  G(\psi)>0$$ by .

It follows from the Slater condition and the theory of Lagrange multipliers [@IK:2008:Lagrange Chapter 1, Theorem 1.6] that there exists a nonnegative number $\lambda$ such that $$\label{eq:VI3}
  b(p,q-p)+(\phi,q-p)+\tau((q(1)-p(1))-\sigma(q(-1)-p(-1))-\lambda\int_I (q-p)dx\geq 0$$ for all $q\in \fK$.

Finally we observe that $$\label{eq:VI4}
 \tilde b(p,q-p)+(F,q-p)+\tau(q(1)-p(1))-\sigma(q(-1)-p(-1))
 \geq 0\qquad\forall\,q\in\fK,$$ where $$\label{eq:taDef}
 \tilde b(q,r)=\int_I q'r'dx+\int_I qr\,dx$$ and $$\label{eq:FDef}
  F=\phi-\lambda-p.$$

The variational inequality defined by , and characterizes the solution of a coercive Neumann obstacle problem on $H^1(I)$. Since $F\in \LT$ and $\psi\in H^2(I)$, we can apply the result in [@Rodrigues:1987:Obstacle Chapter 5, Theorem 3.4] to conclude that $p\in H^2(I)$.

We can deduce the regularity of $(\bar y,\bar u)$ from the relations $p=\bar y'$ and $\bar u=-(\bar y''+f)$.

\[cor:Regularity\] Under the assumption on the data, the solution $(\bar y,\bar u)$ of the optimal control problem – belongs to $H^3(I)\times H^1(I)$.

The result in [@Rodrigues:1987:Obstacle], which is for dimensions $\geq 2$, requires a compatibility condition between $\p\psi/\p n$ and the Neumann boundary condition so that the boundary trace of the normal derivative of the solution of the obstacle problem belongs to the correct Sobolev space. This is not needed in one dimension since the boundary values of the normal derivative are just numbers.

The Karush-Kuhn-Tucker Conditions {#subsec:KKT}
---------------------------------

It follows from , Theorem \[thm:pRegularity\] and integration by parts that $$\label{eq:ObstacleLM}
  b(p,q)+(\phi,q)+\tau q(1)-\sigma q(-1)-\lambda\int_I q\,dx
  +\int_I q\,d\nu=0 \qquad\forall\,q\in H^1(I),$$ where the regular Borel measure $\nu$ is given by $$\label{eq:ObstaclenuDef}
  d\nu=(p''-\phi+\lambda)dx+[p'(-1)+\sigma]d\delta_{-1}-[p'(1)
  +\tau]d\delta_{1},$$ and $\delta_{-1}$ (resp., $\delta_1$) is the Dirac point measure at $-1$ (resp., $1$).

Let $\fA$ be the active set of the derivative constraint , i.e., $$\label{eq:ActiveSet}
  \fA=\{x\in [-1,1]:\,\bar y'(x)=\psi(x)\}=\{x\in[-1,1]:\,p(x)=\psi(x)\}.$$ By a standard argument, $p$ satisfies if and only if $$\label{eq:ObstaclenuProperties}
  \text{$\nu$ is nonnegative and supported on $\fA$.}$$

We can translate – into Karush-Kuhn-Tucker (KKT) conditions for the solution $\bar y'=p\in \cK$ of /, which is summarized in the following theorem.

\[thm:KKT\] There exists a nonnegative number $\lambda$ such that $$\begin{aligned}
 &\int_I p'q'dx+\int_I(\Phi-f')q\,dx+f(1)q(1)-f(-1)q(-1)+
 \int_I q\,d\nu\label{eq:KKT1}\\
 &\hspace{100pt}=\lambda\int_I q\,dx&\qquad\forall\,q\in H^1(I),\notag\\
  &\int_{[-1,1]}(p-\psi)d\nu=0,\label{eq:KKT2}\\[4pt]
  &d\nu=\rho\,dx+\gamma d\delta_{-1}+\zeta d\delta_1,\label{eq:KKT3}\end{aligned}$$ where $$\begin{aligned}
 &\text{$\rho=p''+f'-\Phi+\lambda\in L_2(I)$ is nonnegative a.e.},\label{eq:KKT4}\\
  &\text{$\gamma=p'(-1)+f(-1)$ and $\zeta=-[p'(1)+f(1)]$
  are nonnegative numbers},\label{eq:KKT5}\end{aligned}$$ and $\Phi\in H^1(I)$ satisfies –.

\[rem:KKT\] The KKT conditions – are also sufficient conditions for . Indeed, they imply, for any $q\in\cK$, $$\begin{aligned}
  &\int_I p'(q'-p')dx+\int_I(\Phi-f')q\,dx+f(1)\big(q(1)-p(1)\big)
  -f(-1)\big(q(-1)-p(-1)\big)\\
    &\hspace{50pt}= \lambda\int_I(q-p)dx-\int_I (q-p)d\nu\\
    &\hspace{50pt}=-\int_I (q-\psi)d\nu\geq 0,\end{aligned}$$ which then also implies $\bar y(x)=\int_{-1}^x p(t)dt$ is the solution of .

Finally we observe that Theorem \[thm:KKT\] implies $$\label{eq:KKT}
  \beta\int_I (\bar y''+f)z''dx+\int_I(\bar y-y_d)z\,dx+\int_{[-1,1]}z'd\mu=0
  \qquad\forall\, z\in H^2(I)\cap H^1_0(I),$$ where $$\label{eq:muDef}
\text{$\mu=\beta\nu$ is a nonnegative Borel measure},$$ and $$\label{eq:Complementarity}
  \int_{[-1,1]} (\bar y'-\psi)d\mu=0.$$

The Discrete Problem {#sec:Discrete}
====================

Let $V_h\subset H^2(I)\cap H^1_0(I)$ be the cubic Hermite finite element space (cf. [@Ciarlet:1978:FEM; @BScott:2008:FEM]) associated with a triangulation/partition $\cT_h$ of $I$ with mesh size $h$. The discrete problem is to find $$\label{eq:DiscreteProblem}
  \bar y_h=\argmin_{y_h\in K_h}\frac12\big(\|y_h-y_d\|_\LT^2
  +\beta\|y_h''+f\|_\LT^2\big),$$ where $$\label{eq:KhDef}
  K_h=\{y_h\in V_h:\, P_hy_h'\leq P_h\psi\},$$ and $P_h$ is the interpolation operator associated with the $P_1$ finite element space associated with $\cT_h$, i.e., the constraint is only enforced at the grid points.

The nodal interpolation operator from $C^1([-1,1])$ onto $V_h$ is denoted by $\Pi_h$.

We will use the following standard estimates for $P_h$ and $\Pi_h$ (cf. [@Ciarlet:1978:FEM; @BScott:2008:FEM]) in the error analysis: $$\begin{aligned}
{3}
  \|\zeta- P_h\zeta\|_\LT&\leq Ch|\zeta|_{H^1(I)}
   &\qquad&\forall\,\zeta\in H^1(I),\label{eq:PhError1}\\
   \|\zeta-P_h\zeta\|_\LT&\leq Ch^2|\zeta|_{H^2(I)}
  &\qquad&\forall\,\zeta\in H^2(I),\label{eq:PhError2}\\
 |\zeta-\Pi_h\zeta|_{H^1(I)}+h|\zeta-\Pi_h\zeta|_{H^2(I)}&\leq Ch^2|\zeta|_{H^3(I)}
 &\qquad&\forall\,\zeta\in H^3(I).\label{eq:PihEst}\end{aligned}$$ Here and below we use $C$ to denote a generic positive constant that is independent of the mesh size $h$.

The unique solution $\bar y_h\in K_h$ of the minimization problem defined by and is characterized by the discrete variational inequality $$\beta\int_I (\bar y_h''+f)(y_h''-\bar y_h'')dx
  +\int_I (\bar y_h-y_d)(y_h-\bar y_h)dx\geq0
  \qquad\forall\,y_h\in K_h,$$ which can also be written as $$\label{eq:DVI}
  a(\bar y_h,y_h-\bar y_h)\geq \int_I y_d(y_h-\bar y_h)dx
  -\beta\int_I f(y_h''-\bar y_h'')dx\qquad\forall\,y_h\in K_h,$$ where the bilinear form $a(\cdot,\cdot)$ is defined in .

The error analysis of the finite element method is based on the approach in [@BSung:2017:State] for state constrained optimal control problems that was extended to one dimensional problems with constraints on the derivative of the state in [@BSW:2020:OneD].

We will use the energy norm $\|\cdot\|_a$ defined by $$\label{eq:EnergyNorm}
  \|v\|_a^2=a(v,v)=\|v\|_\LT^2+\beta|v|_{H^2(I)}^2.$$ Note that $$\label{eq:NormEquivalence}
  \|v\|_a\approx \|v\|_{H^2(I)} \qquad\forall\,v\in H^2(I)\cap H^1_0(I)$$ by a Poincaré-Friedrichs inequality [@Necas:2012:Direct].

An Abstract Error Estimate {#subsec:Abstract}
--------------------------

In view of , and the Cauchy-Schwarz inequality, we have $$\begin{aligned}
\label{eq:Preliminary}
  \|\bar y-\bar y_h\|_a^2&=a(\bar y-\bar y_h,\bar y-y_h)+
      a(\bar y-\bar y_h,y_h-\bar y_h)\notag\\
      &\leq \frac12\|\bar y-\bar y_h\|_a^2+\frac12\|\bar y-y_h\|_a^2
      +a(\bar y,y_h-\bar y_h)\\
            &\hspace{60pt}
             -\int_I y_d(y_h-\bar y_h)dx+
     \beta\int_I f(y_h''-\bar y_h'')dx\qquad\forall\,y_h\in K_h.\notag\end{aligned}$$

It follows from , and that $$\begin{aligned}
\label{eq:KKTRelation}
  &a(\bar y,y_h-\bar y_h)-\int_I y_d(y_h-\bar y_h)dx+
     \beta\int_I f(y_h''-\bar y_h'')dx\notag\\
     &\hspace{40pt}=\int_{[-1,1]} (\bar y_h'-y_h')d\mu\notag\\
     &\hspace{40pt}=\int_{[-1,1]} (\bar y_h'-P_h\bar y_h')d\mu+
\int_{[-1,1]} (P_h\bar y_h'-P_h\psi)d\mu+\int_{[-1,1]} (P_h\psi-\psi)d\mu\\
        &\hspace{90pt}+\int_{[-1,1]} (\psi-\bar y')d\mu+
         \int_{[-1,1]} (\bar y'-y_h')d\mu\notag\\
&\hspace{40pt}\leq \int_{[-1,1]} (\bar y_h'-P_h\bar y_h')d\mu
  +\int_{[-1,1]} (P_h\psi-\psi)d\mu
       +\int_{[-1,1]} (\bar y'-y_h')d\mu
       \notag\end{aligned}$$ for all $y_h\in K_h$.

Putting and together, we arrive at the abstract error estimate $$\begin{aligned}
\label{eq:Abstract}
 \|\bar y-\bar y_h\|_a^2&\leq 2\Big(\int_{[-1,1]}
  (\bar y_h'-P_h\bar y_h')d\mu
   +\int_{[-1,1]}(P_h\psi-\psi)d\mu\Big)\\
   &\hspace{50pt} +\inf_{y_h\in K_h}\Big(\|\bar y-y_h\|_a^2+
   2\int_{[-1,1]} (\bar y'-y_h')d\mu\Big).\notag\end{aligned}$$

Concrete Error Estimates {#subsec:Concrete}
------------------------

The three terms on the right-hand side of can be estimated as follows.

First of all, we have $$\begin{aligned}
\label{eq:Est1}
  \int_{[-1,1]}(\bar y_h'-P_h\bar y_h')d\mu&
   =\int_{[-1,1]} \big[(\bar y_h'-\bar y')-P_h(\bar y_h'-\bar y')\big]d\mu+
     \int_{[-1,1]} (\bar y'-P_h\bar y')d\mu\notag\\
 &=\beta \Big(\int_I \big[(\bar y_h'-\bar y')-P_h(\bar y_h'-\bar y')\big]\rho\,dx+
         \int_I (\bar y'-P_h\bar y')\rho\,dx\Big)\\
      &\leq C\big(h\|\bar y-\bar y_h\|_a+h^2|y|_{H^3(I)}\big),\notag\end{aligned}$$ by Corollary \[cor:Regularity\], , , , , and the fact that $\zeta-P_h\zeta$ vanishes at the points $\pm1$ for any $\zeta\in H^1(I)$.

Similarly we can derive $$\label{eq:Est2}
  \int_{[-1,1]}(P_h\psi-\psi)d\mu=\beta\int_I (P_h\psi-\psi)\rho\,dx
  \leq Ch^2|\psi|_{H^2(I)}$$ by and .

Finally we have $$\begin{aligned}
\label{eq:Est3}
  &\inf_{y_h\in K_h}\Big(\|\bar y-y_h\|_a^2+
   2\int_{[-1,1]} (\bar y'-y_h')d\mu\Big)\notag\\
   &\hspace{50pt}\leq \|\bar y-\Pi_h \bar y\|_a^2+2\int_{[-1,1]}
   \big[\bar y'-(\Pi_h\bar y)'\big]d\mu\\
   &\hspace{50pt}= \|\bar y-\Pi_h \bar y\|_a^2+2\beta\int_I
   \big[\bar y'-(\Pi_h\bar y)'\big]\rho\,dx
   \leq Ch^2\big[|\bar y|_{H^3(I)}^2+|\bar y|_{H^3(I)}\big],\notag\end{aligned}$$ by Corollary \[cor:Regularity\], , , , and the fact that $\bar y'-(\Pi_h\bar y)'$ vanishes at $\pm1$.

It follows from – and Young’s inequality that $$\label{eq:EnergyError}
  \|\bar y-\bar y_h\|_a\leq Ch,$$ which immediately implies the following result, where $\bar u_h=-(\bar y_h''+f)$ is the approximation for $\bar u=-(\bar y+f)$.

\[thm:ErrorEstimates\] Under the assumptions on the data in , we have $$|\bar y-\bar y_h|_\HO+\|\bar u-\bar u_h\|_\LT\leq Ch.$$

\[rem:Sharp\] Numerical results in Section \[sec:Numerics\] indicate that the estimate for $\|\bar u-\bar u_h\|_\LT$ in Theorem \[thm:ErrorEstimates\] is sharp.

\[rem:Comparison\] For comparison, the error estimate $$|\bar y-\bar y_h|_\HO+\|\bar u-\bar u_h\|_\LT\leq C_\epsilon h^{\frac12-\epsilon}$$ was obtained in [@BSW:2020:OneD] under the assumptions in .

A Numerical Experiment {#sec:Numerics}
======================

We begin by constructing an example for the problem / with a known exact solution.

An Example {#subsec:Example}
----------

Let $\beta=1$, $$\label{eq:psiExample}
  \psi(x)=\begin{cases}
   1-\frac92 x^2&\qquad -1\leq x\leq 0\\[2pt]
   1&\qquad \phantom{-}0\leq x\leq 1
   \end{cases},$$ and $$\label{eq:ExactSolution}
  \bar y(x)=\int_{-1}^x p(t)dt,$$ where $$\label{eq:pDef}
  p(x)=\begin{cases}
    1-\frac{81}{32}(x-\frac13)^2&\qquad -1\leq x\leq \frac13 \\[2pt]
    1&\qquad\hspace{10pt} \frac13\leq x\leq 1
  \end{cases}.$$

We have $\psi\in H^2(I)$, $$\label{eq:psiExampleConstraint}
  \int_I\psi\, dx =\frac12,$$ $p\in H^2(I)$, $$\label{eq:pSD}
  p''(x)=\begin{cases}
    -81/16& \qquad -1<x<\frac13\\[2pt]
    0&\qquad\hspace{10pt} \frac13<x<1
  \end{cases},$$ $p'(1)=0$, $p'(-1)=27/4$, $$\label{eq:pProperties}
 \int_I p\,dx=0, \quad p\leq\psi \quad\text{and}\quad \fA=\{-1\}\cup[1/3,1].$$

Let $f\in H^1(I)$ be defined by $$\label{eq:fDef}
  f(x)=\begin{cases}\d
    \frac{2}{9\pi}\sin (\pi(3x-1))&\qquad -1< x\leq \frac13 \\[8pt]
    -(x-\frac13)^2 &\qquad \hspace{10pt} \frac13\leq x<1
  \end{cases}.$$ We have $f(-1)=0$, $f(1/3)=0$, $f_-'(1/3)=2/3$, $f'_+(1/3)=0$ and $f(1)=-4/9$. Therefore the function $$\label{eq:PhiExample}
  \Phi(x)=\begin{cases}
    f'(x)&\qquad -1<x< \frac13\\[2pt]
    f'(x)+\frac23&\qquad\hspace{10pt} \frac13< x<1
  \end{cases}$$ belongs to $H^1(I)$ and $$\label{eq:PhiIntegralExample}
  \int_I\Phi\,dx=\int_I f'(x)+\int_\frac13^1\frac23\,dx
  =f(1)-f(-1)+\frac49=0.$$

Finally we take $\lambda=81/16$ and $y_d=\bar y+\Phi'$. Then the KKT conditions – are satisfied with $$\begin{aligned}
  d\nu&=[p''+f'-\Phi+\lambda]dx+[p'(-1)+f(-1)]d\delta_{-1}
  -[p'(1)+f(1)]d\delta_1\\
    &=(211/48)\chi_{[1/3,1]} dx+(27/4)d\delta_{-1}+(4/9)d\delta_1,\end{aligned}$$ where $\chi_{[1/3,1]}$ is the characteristic function of the interval $[1/3,1]$.

Numerical Results {#subsec:Results}
-----------------

We solved the problem in Section \[subsec:Example\] by the finite element method in Section \[sec:Discrete\] on uniform meshes. The results are displayed in Table \[table:Results\].

[|l|c|c|c|c|]{}

------------------------------------------------------------------------

$2/h$ &$\|\bar y-\bar y_h\|_{L_2(I)}$&$\|\bar y-\bar y_h\|_{L_\infty(I)}$ &$|\bar y-\bar y_h|_{H^1(I)}$ &$|\bar y-\bar y_h|_{H^2(I)}$\
&&&&\
$1+2^0$ & 1.430334e-01& 1.625937e-01& 2.581989e-01& 8.660252e-01\
$1+2^1$ & 1.216070e-01& 1.385037e-01& 2.199480e-01& 7.486796e-01\
$1+2^2$ & 4.306657e-02 & 4.679253e-02 & 8.061916e-02 & 4.485156e-01\
$1+2^3$ & 1.613494e-02 & 1.850729e-02 & 2.919318e-02 & 2.573315e-01\
$1+2^4$ & 3.439341e-03 & 3.849954e-03 & 6.315816e-03 & 1.266029e-01\
$1+2^5$ & 9.590453e-04 & 1.087740e-03 & 1.741244e-03 & 6.470514e-02\
$1+2^6$ & 2.256478e-04& 2.542346e-04 &4.125212e-04 & 3.223430e-02\
$1+2^7$ & 5.874304e-05 & 6.639870e-05& 1.067193e-04&1.618687e-02\
$1+2^8$ & 1.425640e-05 &1.608790e-05 &2.549283e-05 & 8.086258e-03\
$1+2^9$ & 3.657433e-06 & 4.124680e-06 & 6.430499e-06 & 4.047165e-03\

We observe $O(h)$ convergence in the $H^2$ norm which agrees with Theorem \[thm:ErrorEstimates\]. On the other hand the convergence in the $H^1$ norm is $O(h^2)$, better than the $O(h)$ convergence predicted by Theorem \[thm:ErrorEstimates\]. The convergence in $L_2$ and $L_\infty$ is also $O(h^2)$.

Concluding Remarks {#sec:Conclusions}
==================

We have shown that higher regularity for the solutions of one dimensional Dirichlet elliptic distributed optimal control problems with pointwise constraints on the derivative of the state can be obtained through a variational inequality satisfied by the derivative of the optimal state. A similar result for one dimensional optimal control problems with mixed boundary conditions was obtained earlier in [@BSW:2020:OneD]. A natural question is: Can these results be extended to higher dimensions?

For analogs of – on a smooth/convex domain $\O\in \R^d$ ($d=2,3$), where $f\in H^1(\O)$ and $\bm{\Psi}\in [H^2(\O)]^d$, one can also derive a variational inequality for the gradient of the optimal state. Observe that the space $\bG$ of the gradients of the states is characterized by (cf. [@GR:1986:NS Chapter I, Section 2.3]) $$\begin{aligned}
 \bG&=\{\nabla y:\,y\in H^2(\O)\cap H^1_0(\O)\}\notag\\
       &=\{\bq\in[ H^1(\O)]^d:\,\Curl\bq=0 \;\text{on}\;\O\;\text{and}
       \;\bn\times\bq=0 \;
       \text{on}\;\p\O\},\end{aligned}$$ where $\bn$ is the unit outward normal along $\p\O$.

Let $\bK$ be the subset of $\bG$ defined by $$\bK=\{\bq\in\bG:\,\bq\leq\bm{\Psi}\;\text{a.e. in $\O$}\}.$$ We assume that $\bK$ is nonempty, which is the case for example if $\bm{\Psi}\geq{\bf 0}$.

The analog of is given by $$\label{eq:DirichletHigher}
  \int_\O  (\bPhi-\nabla f)\cdot(\bq-\bp)dx
  +\int_\O \Div\bp\,\Div(\bq-\bp)\,dx
  +\int_{\p\O} f(\bq-\bp)\cdot\bn\,dS\geq 0$$ for all $\bq\in \bK$, where $\bp=\nabla\bar y\in\bK$, and $\bPhi\in \bG$ is defined by $ \beta\Div\bPhi=y_d-\bar y$, which is an analog of . The variational inequality is uniquely solvable because (cf. [@GR:1986:NS Chapter I, Sections 3.2 and 3.4]) $$\int_\O (\Div \bq)^2dx\geq C_\O|\bq|_{H^1(\O)}^2\qquad \forall\,\bq\in \bG.$$

We can also write as $$\begin{aligned}
\label{eq:DirichletVector}
  &\int_\O  (\bPhi-\nabla f)\cdot(\bq-\bp)dx+
  \int_\O \big[\Div\bp\,\Div(\bq-\bp)
  +\Curl\bp\cdot\Curl(\bq-\bp)\big]\,dx\\
  &\hspace{100pt}+\int_{\p\O} f(\bq-\bp)\cdot\bn\,dS
  \geq 0\qquad\forall\,\bq\in \bK,\notag\end{aligned}$$ which can be interpreted as an obstacle problem for the vector Laplacian operator with natural boundary conditions.

In order to obtain higher regularity for the optimal state $\bar y$, one will need regularity results for /, which unfortunately are not available. Therefore the problem of extending the results in this paper to higher dimensions remains open.

[10]{}

S.C. Brenner, C.B. Davis, and L.-Y. Sung. . , 276:612–626, 2014.

S.C. Brenner, J. Gedicke, and L.-Y. Sung. interior penalty methods for an elliptic distributed optimal control problem on nonconvex polygonal domains with pointwise state constraints. , 56:1758–1785, 2018.

S.C. Brenner, T. Gudi, K. Porwal, and L.-Y. Sung. A [M]{}orley finite element method for an elliptic distributed optimal control problem with pointwise state and control constraints. , 24:1181–1206, 2018.

S.C. Brenner, M. Oh, S. Pollock, K. Porwal, M. Schedensack, and N. Sharma. . In S.C. Brenner, editor, [*Topics in Numerical Partial Differential Equations and Scientific Computing*]{}, volume 160 of [*The IMA Volumes in Mathematics and its Applications*]{}, pages 1–22, Cham-Heidelberg-New York-Dordrecht-London, 2016. Springer.

S.C. Brenner and L.R. Scott. . Springer-Verlag, New York, 2008.

S.C. Brenner and L.-Y. Sung. A new convergence analysis of finite element methods for elliptic distributed optimal control problems with pointwise state constraints. , 55:2289–2304, 2017.

S.C. Brenner, L.-Y. Sung, and Y. Zhang. A quadratic [$C^0$]{} interior penalty method for an elliptic optimal control problem with state constraints. In O. Karakashian X. Feng and Y. Xing, editors, [*Recent Developments in Discontinuous Galerkin Finite Element Methods for Partial Differential Equations*]{}, volume 157 of [*The IMA Volumes in Mathematics and its Applications*]{}, pages 97–132, Cham-Heidelberg-New York-Dordrecht-London, 2013. Springer. (2012 John H. Barrett Memorial Lectures).

S.C. Brenner, L-Y. Sung, and Y. Zhang. interior penalty methods for an elliptic state-constrained optimal control problem with [N]{}eumann boundary condition. , 350:212–232, 2019.

S.C. Brenner, L.Y. Sung, and W. Wollner. Finite element methods for one dimensional elliptic distributed optimal control problems with pointwise constraints on the derivative of the state. , [published online 18 February 2020 (doi.org/10.1007/s11081-020-09491-1)]{}.

E. Casas and J.F. Bonnans. Contr[ô]{}le de syst[è]{}mes elliptiques semilin[é]{}ares comportant des contraintes sur l’[é]{}tat. In H. Brezzis and J.L. Lions, editors, [*Nonlinear Partial Differential Equations and their Applications 8*]{}, pages 69–86. Longman, New York, 1988.

E. Casas and L.A. Fernández. Optimal control of semilinear elliptic equations with pointwise constraints on the gradient of the state. , 27:35–56, 1993.

P.G. Ciarlet. . North-Holland, Amsterdam, 1978.

K. Deckelnick, A. Günther, and M. Hinze. Finite element approximation of elliptic control problems with constraints on the gradient. , 111:335–350, 2009.

I. Ekeland and R. T[é]{}mam. . Classics in Applied Mathematics. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 1999.

V. Girault and P.-A. Raviart. . Springer-Verlag, Berlin, 1986.

W. Gong and N. Yan. A mixed finite element scheme for optimal control problems with pointwise state constraints. , 46:182–203, 2011.

K. Ito and K. Kunisch. . Society for Industrial and Applied Mathematics, Philadelphia, PA, 2008.

D. Kinderlehrer and G. Stampacchia. . Society for Industrial and Applied Mathematics, Philadelphia, 2000.

W. Liu, W. Gong, and N. Yan. A new finite element approximation of a state-constrained optimal control problem. , 27:97–114, 2009.

J. Ne[č]{}as. . Springer, Heidelberg, 2012.

C. Ortner and W. Wollner. A priori error estimates for optimal control problems with pointwise constraints on the gradient of the state. , 118:587–600, 2011.

M. Pierre and J. Sokoł owski. Differentiability of projection and applications. In [*Control of partial differential equations and applications ([L]{}aredo, 1994)*]{}, volume 174 of [*Lecture Notes in Pure and Appl. Math.*]{}, pages 231–240. Dekker, New York, 1996.

J.-F. Rodrigues. . North-Holland Publishing Co., Amsterdam, 1987.

W. Wollner. Optimal control of elliptic equations with pointwise constraints on the gradient of the state in nonsmooth polygonal domains. , 50:2117–2129, 2012.

[^1]: The work of the first and second authors was supported in part by the National Science Foundation under Grant No. DMS-19-13035.
