2pt 4pt

# Douglas–Rachford Algorithm for Control- and State-constrained Optimal
Control Problems

Regina S. Burachik Bethany I. Caldwell11footnotemark: 1 C. Yalçın
Kaya11footnotemark: 1 Mathematics, UniSA STEM, University of South Australia,
Mawson Lakes, S.A. 5095, Australia. Emails<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We consider the application of the Douglas–Rachford (DR) algorithm to solve
linear–quadratic (LQ) control problems with box constraints on the state and
control variables. We split the constraints of the optimal control problem
into two sets: one involving the ODE with boundary conditions, which is
affine, and the other a box. We rewrite the LQ control problems as the
minimization of the sum of two convex functions. We find the proximal mappings
of these functions which we then employ for the projections in the DR
iterations. We propose a numerical algorithm for computing the projection onto
the affine set. We present a conjecture for finding the costates and the state
constraint multipliers of the optimal control problem, which can in turn be
used in verifying the optimality conditions. We carry out numerical
experiments with two constrained optimal control problems to illustrate the
working and the efficiency of the DR algorithm compared to the traditional
approach of direct discretization.

> Key words: Optimal control, Douglas–Rachford algorithm, Harmonic oscillator,
> Control constraints, State constraints, Numerical methods.

Mathematical Subject Classification: 49M37; 49N10; 65K10; 65L10; 34H05

##  1 Introduction

Many of the problems we encounter in the world that can be presented as
optimal control problems contain constraints on both the state and control
variables. Imagine a scheduling problem where the aim is to maximize profits
such as in [26] in which the control variable is the production rate and the
state variable is the inventory level. Naively, one may want to greatly
increase both the inventory level and production rate to maximize profits but
in practice there is a limit to the amount of inventory a factory can hold and
a limit on how quickly humans or machines can work. In order to accurately
model this problem, and many others from a wide rage of application areas such
as manufacturing, engineering, science, economics, etc., we must use an
optimal control problem with both state and control constraints. Although
these problems are commonly found in applications they are much more difficult
to solve than optimal control problems which purely have control constraints.

We focus our attention on linear–quadratic (LQ) state- and control-constrained
optimal control problems. These are infinite dimensional optimization problems
that involve the minimization of a quadratic objective function subject to
linear differential equation (DE) constraints, affine constraints on the state
variables, and affine constraints on the control variables. There is extensive
literature studying LQ control problems as they model many problems from a
wide variety of disciplines – see [14, 16, 1, 27, 28, 15, 24]. Though many of
the LQ control problems posed in the literature contain control constraints
along with the linear DE constraints, it is rarer to see state constraints
included since, as mentioned above, they are much more difficult to deal with.
This paper proposes a unique approach for solving state-constrained problems
by applying the Douglas–Rachford (DR) algorithm.

The DR algorithm is a projection algorithm used to minimize the sum of two
convex functions. In order to apply the algorithm one needs the proximal
operators of the two convex functions. Splitting and projection methods such
as the DR algorithm are a popular area of research in optimization with a
variety of applications – see [20, 2, 9, 3, 6] for their use in sphere packing
problems, protein reconstruction, etc. The use of these methods to solve
discrete-time optimal control problems is not new but there are very few
applications of these methods to continuous-time optimal control problems. The
paper [7] is one in which projection methods were used to solve the energy-
minimizing double integrator problem with control constraints. The papers [11,
12, 13] address more general energy-minimizing LQ control problems with
control constraints. While a collection of general projection methods is used
in [11], the DR algorithm is used in [12, 13].

Following the promising numerical results observed in [12, 13], the present
paper will use the DR algorithm for addressing the more challenging control-
and state-constrained LQ control problems.

The current technique for solving these control- and state-constrained LQ
problems is a direct discretization approach where the infinite-dimensional
problem is reduced to a large scale finite-dimensional problem by use of a
discretization scheme, e.g. a Runge–Kutta method [22, 5]. This discretized
problem is then solved through the use of an optimization modelling language
such as AMPL [19] paired with a (finite-dimensional) large-scale optimization
software such as Ipopt [33].

In the present paper (as was done in [7, 11] for control constraints only) we
make the following contributions for the numerical solution of LQ control
problems with state and control constraints.

  * (i)

We rewrite the original LQ control problem as the minimization of the sum of
two convex functions.

  * (ii)

We find the proximal mappings of these functions by solving the resulting
calculus of variations subproblems (Theorems 1–2).

  * (iii)

We define and apply the DR algorithm to the constrained LQ control problem
(Algorithm 1).

  * (iv)

We propose an algorithm for computing the projection onto the affine
constraint set defined by the ODEs (Algorithm 2).

  * (v)

We carry out numerical experiments to illustrate the efficiency of the DR
algorithm over the traditional approach of direct discretization.

In addition to these, we also make the following contribution, which is new
even for the purely control-constrained LQ control problems.

  * (vi)

We present a conjecture for finding the costates. Moreover, we present a
technique for finding the multiplier for the state constraints in the case
only one state constraint is active at a given time. We effectively employ the
costates and the state control multiplier for a numerical verification of the
optimality conditions.

The two convex functions mentioned above are defined as follows: the
information from the ODE constraints appears in one function while the
information from the control and state constraints along with the objective
functional appear in the other. We derive the proximal mappings of the two
convex functions without reducing the original infinite-dimensional
optimization problem to a finite-dimensional one, though we need to discretize
the state and control variables over a partition of time when implementing the
DR algorithm since a digital computer cannot iterate with functions.

The paper is structured as follows. In Section 2 we give the problem
formulation and optimality conditions for the optimal control problem. In
Section 3 we derive the proximal mappings used in the implementation of the DR
algorithm. Section 4 introduces the DR algorithm. Then Section 5 begins with
an algorithm for computing one of the proximal mappings (namely the projection
onto an affine set) and a conjecture used to obtain the costate variable of
the original LQ control problem. Next, this section introduces two example
problems, a harmonic oscillator and a mass–spring system. At the end of this
section, numerical experiments for the DR algorithm and AMPL–Ipopt suite and
their comparisons are given for these two problems. Finally, concluding
remarks and comments for future work are provided in Section 6.

##  2 Optimal Control Problem

In this section we formulate the general optimal control problem that will be
the focus of this paper. We give some necessary definitions and provide
conditions for optimality from optimal control theory.

Before introducing the optimal control problem we will give some standard
definitions. Unless otherwise stated all vectors are column vectors. Let
${\cal L}^{2}([t_{0},t_{f}];\mathbb{R}^{q})$ be the Banach space of Lebesgue
measurable functions $z:[t_{0},t_{f}]\rightarrow\mathbb{R}^{q}$, with finite
${\cal L}^{2}$ norm, namely,

${\cal
L}^{2}([t_{0},t_{f}];\mathbb{R}^{q}):=\left\\{z:[t_{0},t_{f}]\rightarrow\mathbb{R}^{q}\,\>|\,\>\|z\|_{{\cal
L}^{2}}:=\left(\int_{t_{0}}^{t_{f}}\|z(t)\|^{2}\,dt\right)^{1/2}<\infty\right\\}$

where $\|\cdot\|$ is the $\ell_{2}$ norm in $\mathbb{R}^{q}$. Furthermore,
${\cal W}^{1,2}([t_{0},t_{f}];\mathbb{R}^{q})$ is the Sobolev space of
absolutely continuous functions, namely

${\cal W}^{1,2}([t_{0},t_{f}];\mathbb{R}^{q}):=\left\\{z\in{\cal
L}^{2}([t_{0},t_{f}];\mathbb{R}^{q})\,|\,\dot{z}:=dz/dt\in{\cal
L}^{2}([t_{0},t_{f}];\mathbb{R}^{q})\right\\},$

endowed with the norm

$\|z\|_{{\cal W}^{1,2}}:=(\|z\|_{{\cal L}^{2}}^{2}+\|\dot{z}\|_{{\cal
L}^{2}}^{2})^{1/2}.$

With these definitions we can state now the general LQ control problem as
follows.

$\displaystyle\mbox{(P) }\left\\{\begin{array}[]{rl}\displaystyle\min&\ \
\displaystyle\frac{1}{2}\int_{t_{0}}^{t_{f}}\left(x(t)^{T}Q(t)x(t)+u(t)^{T}R(t)u(t)\right)\,dt\\\\[14.22636pt]
\mbox{subject to}&\ \ \dot{x}(t)=A(t)x(t)+B(t)u(t)\,,\ \ x(t_{0})=x_{0}\,,\ \
x(t_{f})=x_{f}\,,\\\\[5.69054pt] &\ \ \underline{u}\leq
u(t)\leq\overline{u}\,,\ \ \underline{x}\leq
x(t)\leq\overline{x}\,.\end{array}\right.$

The state variable $x\in{\cal W}^{1,2}([t_{0},t_{f}];\mathbb{R}^{n})$, with
$x(t):=(x_{1}(t),\ldots,x_{n}(t))\in\mathbb{R}^{n}$, and the control variable
$u\in{\cal L}^{2}([t_{0},t_{f}];\mathbb{R}^{m})$, with
$u(t):=(u_{1}(t),\ldots,u_{m}(t))\in\mathbb{R}^{m}$. The time varying matrices
$A:[t_{0},t_{f}]\rightarrow\mathbb{R}^{n\times n}$,
$B:[t_{0},t_{f}]\to\mathbb{R}^{n\times m}$,
$Q:[t_{0},t_{f}]\rightarrow\mathbb{R}^{n\times n}$ and
$R:[t_{0},t_{f}]\rightarrow\mathbb{R}^{m\times m}$ are continuous. For every
$t\in[t_{0},t_{f}]$, the matrices $Q(t)$ and $R(t)$ are symmetric, and
respectively positive semi-definite and positive definite. For clarity of
arguments, these matrices are assumed to be diagonal, namely
$Q(t):=\operatorname*{diag}(q_{1}(t),\ldots,q_{n}(t))$ and
$R(t):=\operatorname*{diag}(r_{1}(t),\ldots,r_{m}(t))$. These diagonality
assumptions particularly simplify proximal mapping expressions that appear
later. The initial and terminal states are given as $x_{0}$ and $x_{f}$
respectively.

###  2.1 Optimality conditions

In this section we state the Maximum Principle by using the direct adjoining
approach from [21]. We start by defining the extended Hamiltonian function
$H:\mathbb{R}^{n}\times\mathbb{R}^{m}\times\mathbb{R}^{n}\times\mathbb{R}^{n}\times\mathbb{R}^{n}\times[t_{0},t_{f}]\to\mathbb{R}$
for Problem (P) as

$H(x(t),u(t),\lambda(t),\mu^{1}(t),\mu^{2}(t),t):=\frac{1}{2}(x(t)^{T}Q(t)x(t)+u(t)^{T}R(t)u(t))+\lambda(t)^{T}(A(t)x(t)\\\
+B(t)u(t))+\mu^{1}(t)^{T}(x(t)-\overline{x}(t))+\mu^{2}(t)^{T}(\underline{x}(t)-x(t)),$

where the adjoint variable vector
$\lambda:[t_{0},t_{f}]\rightarrow\mathbb{R}^{n}$ with
$\lambda(t):=(\lambda_{1}(t),\dots,\lambda_{n}(t))\in\mathbb{R}^{n}$ and the
state constraint multiplier vectors
$\mu^{1},\mu^{2}:[t_{0},t_{f}]\rightarrow\mathbb{R}^{n}$ with
$\mu^{i}(t):=(\mu_{1}^{i}(t),\dots,\mu_{n}^{i}(t))\in\mathbb{R}^{n}$, $i=1,2$.
For brevity, we use the following notation,

$H[t]:=H(x(t),u(t),\lambda(t),\mu^{1}(t),\mu^{2}(t),t)\,.$

The adjoint variable vector is assumed to satisfy

$\dot{\lambda}(t):=-H_{x}[t]=-Q(t)x(t)-A(t)^{T}\lambda(t)-\mu^{1}(t)+\mu^{2}(t)$
(1)

for a.e. $t\in[t_{0},t_{f}]$, where $H_{x}:=\partial H/\partial x$.

Maximum Principle. Suppose the pair

$(x,u)\in{\cal W}^{1,2}([t_{0},t_{f}];\mathbb{R}^{n})\times{\cal
L}^{2}([t_{0},t_{f}];\mathbb{R}^{m})$

is optimal for Problem (P). Then there exists a piecewise continuous adjoint
variable vector $\lambda\in{\cal W}^{1,2}([t_{0},t_{f}];\mathbb{R}^{n})$ as
defined in (1) and piecewise continuous multipliers $\mu^{1},\mu^{2}\in{\cal
L}^{2}([t_{0},t_{f}];\mathbb{R}^{n})$, such that
$(\lambda(t),\mu^{1}(t),\mu^{2}(t))\neq\mathbf{0}$ for all
$t\in[t_{0},t_{f}]$, and that for a.e. $t\in[t_{0},t_{f}]$,

$\displaystyle u_{i}(t)$
$\displaystyle=\operatorname*{argmin}_{\underline{u}_{i}\leq
v_{i}\leq\overline{u}_{i}}H(x(t),u_{1}(t),\ldots,v_{i},\dots,u_{m}(t),\lambda(t),\mu^{1}(t),\mu^{2}(t),t)$
$\displaystyle=\operatorname*{argmin}_{\underline{u}_{i}\leq
v_{i}\leq\overline{u}_{i}}\frac{1}{2}r_{i}(t)v_{i}^{2}+\lambda(t)^{T}b_{i}(t)v_{i}\,,$
(2)

for $i=1,\ldots,m$, where $b_{i}(t)$ is the $i$th column of the matrix $B(t)$
and $r_{i}(t)$ is the $i$th diagonal element of $R(t)$. Moreover, the
multipliers $\mu^{1}(t),\mu^{2}(t)$ must satisfy the complementarity
conditions

$\displaystyle\mu_{i}^{1}(t)\geq
0,\qquad\mu_{i}^{1}(t)(x_{i}(t)-\overline{x}_{i}(t))=0$ (3)
$\displaystyle\mu_{i}^{2}(t)\geq
0,\qquad\mu_{i}^{2}(t)(\underline{x}_{i}(t)-x_{i}(t))=0$ (4)

for all $i=1,\dots,n$.

Suppose $\underline{u}_{i}=-\infty$, $\overline{u}_{i}=\infty$,
$i=1,\ldots,m$, i.e., the control vector is unconstrained. Then (2) becomes

$H_{u_{i}}[t]=0,$

i.e.,

$r_{i}(t)u_{i}(t)+b_{i}(t)^{T}\lambda(t)=0\,,$ (5)

$i=1,\ldots,m$. Then (5) can be solved for $u_{i}(t)$ as

$u_{i}(t)=-\frac{1}{r_{i}(t)}b_{i}(t)^{T}\lambda(t)\,,$ (6)

for $i=1,\ldots,m$; or using the matrices $B(t),R(t)$,

$u(t)=-[R(t)]^{-1}B(t)^{T}\lambda(t)\,.$ (7)

When we consider the constraints on $u$, one gets from (2),

$\displaystyle u_{j}(t)=\left\\{\begin{array}[]{rl}\underline{u}_{j},&\mbox{if
}\frac{1}{r_{j}(t)}b_{j}(t)^{T}\lambda(t)\geq\underline{u}_{j},\\\\[5.69054pt]
-\frac{1}{r_{j}(t)}b_{j}(t)^{T}\lambda(t),&\mbox{if
}\underline{u}_{j}\leq-\frac{1}{r_{j}(t)}b_{j}(t)^{T}\lambda(t)\leq\overline{u}_{j},\\\\[5.69054pt]
\overline{u}_{j},&\mbox{if
}\frac{1}{r_{j}(t)}b_{j}(t)^{T}\lambda(t)\leq\overline{u}_{j},\end{array}\right.$
(11)

for a.e. $t\in[t_{0},t_{f}]$, $j=1,\dots,m$.

##  3 Splitting and Proximal Mappings

In this section, we rewrite Problem (P) as the minimization of the sum of two
convex functions $f$ and $g$, and give the proximal mappings for these
functions in Theorems 1–2.

We split the constraints from (P) into two sets $\cal{A},\cal{B}$ given as

$\displaystyle{\cal A}:=$ $\displaystyle\big{\\{}(x,u)\in{\cal
L}^{2}([t_{0},t_{f}];\mathbb{R}^{n})\times{\cal
L}^{2}([t_{0},t_{f}];\mathbb{R}^{m})\ |\mbox{ which solves }$ (12)
$\displaystyle\dot{x}(t)=A(t)x(t)+B(t)u(t),\ \ x(t_{0})=x_{0},\ \
x(t_{f})=x_{f},\mbox{ for all }t\in[t_{0},t_{f}]\big{\\}}\,,$
$\displaystyle{\cal B}:=$ $\displaystyle\big{\\{}(x,u)\in{\cal
L}^{2}([t_{0},t_{f}];\mathbb{R}^{n})\times{\cal
L}^{2}([t_{0},t_{f}];\mathbb{R}^{m})\ |\ \underline{x}_{i}\leq
x_{i}(t)\leq\overline{x}_{i}\,,\ \underline{u}_{j}\leq
u_{j}(t)\leq\overline{u}_{j}\,,$ (13) $\displaystyle\mbox{for all
}t\in[t_{0},t_{f}],\ i=1,\dots,n,\ j=1,\dots,m\big{\\}}\,.$

Despite previously defining $x\in{\cal W}^{1,2}([t_{0},t_{f}];\mathbb{R}^{n})$
in our sets ${\cal A},{\cal B}$ we let $x\in{\cal
L}^{2}([t_{0},t_{f}];\mathbb{R}^{n})$ to simplify the calculation of the
proximal mappings. By definition of ${\cal
W}^{1,2}([t_{0},t_{f}];\mathbb{R}^{n})$ if $x\in{\cal A}$ then $x\in{\cal
W}^{1,2}([t_{0},t_{f}];\mathbb{R}^{n})$.

We assume that the control system $\dot{x}(t)=A(t)x(t)+B(t)u(t)$ is
controllable; in other words the control system can be driven from any $x_{0}$
to any other $x_{f}$—for a precise definition of controllability and the tests
for controllability, see [29]. Then there exists a (possibly not unique)
$u(\cdot)$ such that, when this $u(\cdot)$ is substituted, the boundary-value
problem given in the expression for ${\cal A}$ has a solution $x(\cdot)$. In
other words, ${\cal A}\neq\emptyset$. Also, clearly, ${\cal B}\neq\emptyset$.
We note that the constraint set ${\cal A}$ is an affine subspace. Given that
${\cal B}$ is a box, the constraints turn out to be two convex sets in Hilbert
space. Since every sequence converging in ${\cal L}^{2}$ has a subsequence
converging pointwise, it is straightforward to see that the set ${\cal B}$ is
closed in ${\cal L}^{2}$. The closedness of ${\cal A}$ will be established
later on as a consequence of Theorem 2 (see Remark 2).

Fix $\beta>0$ and let

$f(x,u):=\iota_{\cal{B}}(x,u)+\frac{\beta}{2}(x(t)^{T}Q(t)x(t)+u(t)^{T}R(t)u(t))\quad\mbox{and}\quad
g(x,u):=\iota_{\cal{A}}(x,u)\,,$ (14)

where $\iota_{\cal{C}}$ is the indicator function of the set $\cal{C}$, namely

$\iota_{\cal C}(x,u):=\left\\{\begin{array}[]{ll}0\,,&\mbox{ if\ \
}(x,u)\in{\cal C}\,,\\\ \infty\,,&\mbox{ otherwise}.\end{array}\right.$

Problem (P) is then equivalent to the following problem.

$\min\ \ f(x,u)+g(x,u).$ (15)

In our setting, we assume that we are able to compute the projector operators
$P_{\cal{A}}$ and $P_{\cal{B}}$. These operators project a given point onto
each of the constraint sets $\cal{A}$ and $\cal{B}$, respectively. Recall that
the proximal mapping of a functional $h$ is defined by [8, Definition 24.1].
For our setting,

$\operatorname{Prox}_{h}(x,u):=\operatorname*{argmin}_{y\in{\cal
L}^{2}([t_{0},t_{f}];\mathbb{R}^{n})\atop v\in{\cal
L}^{2}([t_{0},t_{f}];\mathbb{R}^{m})}\left(h(y,v)+\frac{1}{2}\|y-x\|_{{\cal
L}^{2}}^{2}+\frac{1}{2}\|v-u\|_{{\cal L}^{2}}^{2}\right),$ (16)

for any $(x,u)\in{\cal L}^{2}([t_{0},t_{f}];\mathbb{R}^{n})\times{\cal
L}^{2}([t_{0},t_{f}];\mathbb{R}^{m})$. Recall that the projection $P_{\cal
C}(x,u)$ of a point $(x,u)$ onto ${\cal C}$ is characterized by $P_{\cal
C}(x,u)\in{\cal C}$ and, $\forall(y,v)\in{\cal C}$, $\langle(y,v)-P_{\cal
C}(x,u)|(x,u)-P_{\cal C}(x,u)\rangle\leq 0$ [8, Theorem 3.16].

Note that $\operatorname{Prox}_{\iota_{\cal C}}=P_{\cal C}$.

In order to implement the Douglas–Rachford algorithm we must write the
proximal mappings $f$ and $g$. The proofs of Theorems 1 and 2 below follow the
broad lines of proof of Lemma 2 in [13]. In both theorems, the major
difference from [13] is that the proximal operators in the current paper have
two variables $x^{-}$ and $u^{-}$. Thanks to separability, the proof of
Theorem 1 is a straightforward modification of the corresponding part of the
proof of [13, Lemma 2]. We include a full proof of Theorem 1 for the
convenience of the reader. On the other hand, the proof of Theorem 2 deals
with the solution of a more involved optimal control subproblem, namely
Problem (Pg).

######  Theorem 1.

The proximal mapping of $f$ is given as
$\operatorname{Prox}_{f}(x^{-},u^{-})=(y,v)$ such that the components of $y$
and $v$ are expressed as

$\displaystyle y_{i}(t)=\left\\{\begin{array}[]{rl}\overline{x}_{i},&\mbox{if
}\frac{1}{\beta q_{i}+1}x^{-}_{i}(t)\geq\overline{x}_{i},\\\\[5.69054pt]
\frac{1}{\beta q_{i}+1}x^{-}_{i}(t),&\mbox{if
}\underline{x}_{i}\leq\frac{1}{\beta
q_{i}+1}x^{-}_{i}(t)\leq\overline{x}_{i},\\\\[5.69054pt]
\underline{x}_{i},&\mbox{if }\frac{1}{\beta
q_{i}+1}x^{-}_{i}(t)\leq\underline{x}_{i},\end{array}\right.$ (20)
$\displaystyle v_{j}(t)=\left\\{\begin{array}[]{rl}\overline{u}_{j},&\mbox{if
}\frac{1}{\beta r_{j}+1}u^{-}_{j}(t)\geq\overline{u}_{j},\\\\[5.69054pt]
\frac{1}{\beta r_{j}+1}u^{-}_{j}(t),&\mbox{if
}\underline{u}_{j}\leq\frac{1}{\beta
r_{j}+1}u^{-}_{j}(t)\leq\overline{u}_{j},\\\\[5.69054pt]
\underline{u}_{j},&\mbox{if }\frac{1}{\beta
r_{j}+1}u^{-}_{j}(t)\leq\underline{u}_{j},\end{array}\right.$ (24)

for all $t\in[t_{0},t_{f}],\ i=1,\dots,n,\ j=1,\dots,m$.

###### Proof.

From (16) and the definition of $f$ in (14) we have that

$\displaystyle\operatorname{Prox}_{f}(x^{-},u^{-})=$
$\displaystyle\operatorname*{argmin}_{x,u}\iota_{\cal{B}}(x,u)+\frac{1}{2}\int_{t_{0}}^{t_{f}}\Bigl{(}\beta
x(t)^{T}Q(t)x(t)+\|x(t)-x^{-}(t)\|^{2}$ $\displaystyle+\beta
u(t)^{T}R(t)u(t)+\|u(t)-u^{-}(t)\|^{2}\Bigr{)}\,dt.$

In other words, to find $\operatorname{Prox}_{f}(x^{-},u^{-})$ we need to find
$(y,v)$ that solves

$\displaystyle\mbox{(Pf) }\left\\{\begin{array}[]{rl}\displaystyle\min&\ \
\displaystyle\frac{1}{2}\int_{t_{0}}^{t_{f}}\Bigl{(}\beta
y(t)^{T}Q(t)y(t)+\|y(t)-x^{-}(t)\|^{2}\\\\[11.38109pt] &\hskip
28.45274pt+\beta
v(t)^{T}R(t)v(t)+\|v(t)-u^{-}(t)\|^{2}\Bigr{)}\,dt\\\\[14.22636pt]
\mbox{subject to}&\ \ \underline{x}_{i}\leq y_{i}(t)\leq\overline{x}_{i}\,,\
\underline{u}_{j}\leq v_{j}(t)\leq\overline{u}_{j}\\\\[5.69054pt] &\ \
\mbox{for all $t\in[t_{0},t_{f}],\ i=1,\dots,n,\
j=1,\dots,m$.}\end{array}\right.$

Problem (Pf) is separable in the variables $y$ and $v$ so we can consider the
problems of minimizing w.r.t. $y$ and $v$ individually and thus solve the two
subproblems

$\displaystyle\mbox{(Pf1) }\left\\{\begin{array}[]{rl}\displaystyle\min&\ \
\displaystyle\frac{1}{2}\int_{t_{0}}^{t_{f}}\beta
y(t)^{T}Q(t)y(t)+\|y(t)-x^{-}(t)\|^{2}\,dt\\\\[14.22636pt] \mbox{subject to}&\
\ \underline{x}_{i}\leq y_{i}(t)\leq\overline{x}_{i}\,,\ \ \mbox{for all
$t\in[t_{0},t_{f}],\ i=1,\dots,n$,}\end{array}\right.$

and

$\displaystyle\mbox{(Pf2) }\left\\{\begin{array}[]{rl}\displaystyle\min&\ \
\displaystyle\frac{1}{2}\int_{t_{0}}^{t_{f}}\beta
v(t)^{T}R(t)v(t)+\|v(t)-u^{-}(t)\|^{2}\,dt\\\\[14.22636pt] \mbox{subject to}&\
\ \underline{u}_{j}\leq v_{j}(t)\leq\overline{u}_{j}\ \ \mbox{for all
$t\in[t_{0},t_{f}],\ j=1,\dots,m$.}\end{array}\right.$

The solution to Problem (Pf1) is given by

$y_{i}(t)=\operatorname*{argmin}_{\underline{x}_{i}\leq
z_{i}\leq\overline{x}_{i}}\big{(}\,\beta\,q_{i}\,z_{i}^{2}+(z_{i}-x_{i}^{-}(t))^{2}\,\big{)}\,,$

$i=1,\ldots,n$, which, after straightforward manipulations, yields (20). The
solution to Problem (Pf2) is obtained similarly as

$v_{j}(t)=\operatorname*{argmin}_{\underline{u}_{j}\leq
w_{j}\leq\overline{u}_{j}}\big{(}\,\beta\,r_{j}\,w_{j}^{2}+(w_{j}-u_{j}^{-}(t))^{2}\,\big{)}\,,$

$j=1,\ldots,m$, which yields (24) after straightforward manipulations. ∎

######  Theorem 2.

The proximal mapping of $g$ is given as
$\operatorname{Prox}_{g}(x^{-},u^{-})=P_{\cal{A}}(x^{-},u^{-})=(y,v)$ such
that

$\displaystyle y(t)=$ $\displaystyle\ x(t),$ (25) $\displaystyle v(t)=$
$\displaystyle\ u^{-}(t)-B(t)^{T}\lambda(t),$ (26)

where $x(t),\lambda(t)$ are obtained by solving the two-point boundary-value
problem (TPBVP)

$\begin{array}[]{ll}&\dot{x}(t)=A(t)x(t)+B(t)u^{-}(t)-B(t)B(t)^{T}\lambda(t)\,,\
\ x(t_{0})=x_{0}\,,\ \ x(t_{f})=x_{f}\,,\\\\[5.69054pt]
&\dot{\lambda}(t)=-x(t)+x^{-}(t)-A(t)^{T}\lambda(t)\,.\end{array}$ (27)

###### Proof.

Using [8, Example 12.25], and the definition of $g$ in (14),

$\operatorname{Prox}_{g}(x^{-},u^{-})=\operatorname{Prox}_{\iota_{\cal{A}}}(x^{-},u^{-})=P_{\cal{A}}(x^{-},u^{-}),$

which verifies the very first assertion. From (16) and the definition of $g$
in (14) we have that

$\operatorname{Prox}_{g}(x^{-},u^{-})=\operatorname*{argmin}_{x,u}\iota_{\cal{A}}(x,u)+\frac{1}{2}\int_{t_{0}}^{t_{f}}(\|x(t)-x^{-}(t)\|^{2}+\|u(t)-u^{-}(t)\|^{2})\,dt.$

In other words, to find $\operatorname{Prox}_{g}(x^{-},u^{-})$ we need to find
$(y,v)$ that solves the problem

$\displaystyle\mbox{(Pg) }\left\\{\begin{array}[]{rl}\displaystyle\min&\ \
\displaystyle\frac{1}{2}\int_{t_{0}}^{t_{f}}\|y(t)-x^{-}(t)\|^{2}+\|v(t)-u^{-}(t)\|^{2}\,dt\\\\[14.22636pt]
\mbox{subject to}&\ \ \dot{y}(t)=A(t)y(t)+B(t)v(t)\,,\ \ y(t_{0})=x_{0}\,,\ \
y(t_{f})=x_{f}\,.\end{array}\right.$

Problem (Pg) is an optimal control problem where $y(t)$ is the state variable
and $v(t)$ is the control variable. The Hamiltonian for Problem (Pg) is

$H(y(t),v(t),\lambda(t),t):=\frac{1}{2}(\|y(t)-x^{-}(t)\|^{2}+\|v(t)-u^{-}(t)\|^{2})+\lambda(t)^{T}(A(t)y(t)+B(t)v(t))$

and the associated costate equation is

$\dot{\lambda}(t)=-\partial H/\partial y=-y(t)+x^{-}(t)-A(t)^{T}\lambda(t).$
(28)

If $v$ is the optimal control for Problem (Pg) then, by the maximum principle,
$H_{v}[t]=0$ for all $t\in[t_{0},t_{f}]$. In other words,

$v(t)-u^{-}(t)+B(t)^{T}\lambda(t)=0\,,$

for all $t\in[t_{0},t_{f}]$, a re-arrangement of which yields (26). Collecting
together the ODE in Problem (Pg) and the ODE in (28), substituting $v(t)$ from
(26), and assigning $y(t)=x(t)$, result in the TPBVP in (27). ∎

######  Remark 1.

We note from Theorem 2 that $\operatorname{Prox}_{g}$ is the projection onto
the affine set $\cal A$. Unlike $\operatorname{Prox}_{f}$, in general, we
cannot find an analytical solution to (27) to obtain $\operatorname{Prox}_{g}$
(or the projection onto $\cal A$) so for this purpose we will propose
Algorithm 2 in Section 4.1 to get a numerical solution. $\Box$

######  Remark 2.

From Theorem 2 we see that every pair $(x,u)\in{\cal L}^{2}\times{\cal L}^{2}$
has a projection onto ${\cal A}$. In other words, ${\cal A}$ is a Chebyshev
set. It is well-known that every Chebyshev set is closed (see [8, Remark
3.11(i)]). Hence, the set ${\cal A}$ is closed in the topology of ${\cal
L}^{2}\times{\cal L}^{2}$. $\Box$

##  4 Douglas–Rachford Algorithm

The application of the Douglas–Rachford (DR) algorithm to our problem is
slightly different to that in [7, 12]. Since we are solving control- and
state-constrained optimal control problems we must define the proximal
mappings at the pair $(x,u)$, rather than just at $u$ as in [7, 12]. Thus in
the implementation of the DR algorithm we give the iterations for the pair of
iterates $(x^{k},u^{k})$ rather than for $u^{k}$ alone.

Given $\beta>0$, we specialize the DR algorithm (see [17], [25] and [18]) to
the case of minimizing the sum of the two functions $f$ and $g$ as in
(14)–(15). The DR operator associated with the ordered pair $(f,g)$ is defined
by

$T:=\operatorname{Id}-\operatorname{Prox}_{f}+\operatorname{Prox}_{g}(2\operatorname{Prox}_{f}-\operatorname{Id})\,.$

Application of the operator to our case is given by

$\displaystyle T(x,u)$
$\displaystyle=(x,u)-\operatorname{Prox}_{f}(x,u)+P_{\cal{A}}(2\operatorname{Prox}_{f}(x,u)-(x,u))\,,$
(29)

where the proximal mappings of $f$ and $g$ are provided as in Theorems 1–2.
Let $X$ be an arbitrary Hilbert space. Now fix $x_{0}\in X$. Given $x_{n}\in
X$, $k\geq 0$, the DR iterations are set as follows.

$\displaystyle(b_{x,k},b_{u,k})$
$\displaystyle:=\operatorname{Prox}_{f}(x_{k},u_{k})\,,$
$\displaystyle(x_{k+1},u_{k+1})$
$\displaystyle:=T(x_{k},u_{k})=(x_{k},u_{k})-(b_{x,k},b_{u,k})+P_{\cal{A}}(2(b_{x,k},b_{u,k})-(x_{k},u_{k}))\,.$

The DR algorithm is implemented as follows. We define a new parameter
$\gamma:=1/(1+\beta)$ where $\beta$ is the parameter multiplying the objective
as in (14) and Theorem 1. The choice of $\gamma\in]0,1[$ can be made because
changing $\beta$ does not affect the solution of Problem (P).

######  Algorithm 1.

(Douglas–Rachford)

Step 1

(Initialization) Choose a parameter $\gamma\in\left]0,1\right[$ and the
initial iterate $(x^{0},u^{0})$ arbitrarily. Choose a small parameter
$\varepsilon>0$, and set $k=0$.

Step 2

(Projection onto ${\cal B}$) Set $(x^{-},u^{-})=\gamma\,(x^{k},u^{k})$.
Compute $(\widetilde{x},\widetilde{u})=P_{{\cal B}}(x^{-},u^{-})$ by using
(20)–(24).

Step 3

(Projection onto ${\cal A}$) Set
$(x^{-},u^{-}):=2\,(\widetilde{x},\widetilde{u})-(x^{k},u^{k})$. Compute
$(\widehat{x},\widehat{u})=P_{{\cal A}}(x^{-},u^{-})$ by using (25)–(26) or
Algorithm 2.

Step 4

(Update) Set
$(x^{k+1},u^{k+1}):=(x^{k},u^{k})+(\widehat{x},\widehat{u})-(\widetilde{x},\widetilde{u})$.

Step 5

(Stopping criterion) If
$\max(\|x^{k+1}-x^{k}\|_{L^{\infty}},\|u^{k+1}-u^{k}\|_{L^{\infty}})\leq\varepsilon$
or $k>=200$, then return $(\widetilde{x},\widetilde{u})$ and stop. Otherwise,
set $k:=k+1$ and go to Step 2.

######  Remark 3.

We point out that, in general, only weak convergence is guaranteed for the DR
algorithm (see [32, Theorem 1] or [10, Theorem 4.4]). It should be noted that
the final iterate of the state–control pair in the box $\cal B$ is returned by
the algorithm. $\Box$

###  4.1 Algorithm for projector onto $\cal A$

We introduce a procedure for numerically projecting onto $\cal A$ that is an
extension of Algorithm 2 from [12] to the case of LQ control problems with
state and control constraints. The procedure below (Algorithm 2) can be
employed in Step 3 of the DR algorithm. In the procedure we effectively solve
the TPBVP in (27) by implementing the standard shooting method [4, 23, 31].
Throughout the steps of Algorithm 2 below, we will solve the ODEs in (27),
rewritten here in matrix form as

$\displaystyle\begin{bmatrix}\dot{x}(t)\\\ \dot{\lambda}^{DR}(t)\end{bmatrix}$
$\displaystyle=\begin{bmatrix}A(t)&-B(t)B(t)^{T}\\\ -I_{n\times
n}&-A(t)^{T}\end{bmatrix}\begin{bmatrix}x(t)\\\
\lambda^{DR}(t)\end{bmatrix}+\begin{bmatrix}0_{n\times n}&B(t)\\\ I_{n\times
n}&0_{n\times m}\end{bmatrix}\begin{bmatrix}x^{-}(t)\\\
u^{-}(t)\end{bmatrix},$ (30)

with various initial conditions (IC):

$\displaystyle\mbox{(i)}\begin{bmatrix}x(t_{0})\\\
\lambda^{DR}(t_{0})\end{bmatrix}=\begin{bmatrix}x_{0}\\\ 0\end{bmatrix},\ \
\mbox{(ii)}\begin{bmatrix}x(t_{0})\\\
\lambda^{DR}(t_{0})\end{bmatrix}=\begin{bmatrix}x_{0}\\\ e_{i}\end{bmatrix},\
\ \mbox{(iii)}\begin{bmatrix}x(t_{0})\\\
\lambda^{DR}(t_{0})\end{bmatrix}=\begin{bmatrix}x_{0}\\\
\lambda^{DR}_{0}\end{bmatrix}.$ (31)

In the above equations, we are using $\lambda^{DR}$, instead of just
$\lambda$, to emphasize the fact that $\lambda^{DR}$ is the costate variable
emanating from solving Problem (Pg) to compute the projection onto $\cal A$
within the DR algorithm. We reiterate that Problem (Pg) is more involved than
its counterpart in [13], which leads to the ODE in (30) which in turn is more
complicated than its counterpart in [12].

######  Algorithm 2.

(Numerical Computation of the Projector onto ${\cal A}$)

Step 0

(Initialization) The following are given: Current iterate $u^{-}$, the system
and control matrices $A(t)$ and $B(t)$, the numbers of state and control
variables $n$ and $m$, and the initial and terminal states $x_{0}$ and
$x_{f}$, respectively. Define $z(t,\lambda_{0}):=x(t)$.

Step 1

(Near-miss function) Solve (30) with ICs in (31)(i) to find
$z(t_{f},0)=x(t_{f})$.
Set $\varphi(0):=z(t_{f},0)-x_{f}$.

Step 2

(Jacobian) For $i=1,\ldots,n$, solve (30) with ICs in (31)(ii), to get
$z(t_{f},e_{i})$.
Set $\beta_{i}(t):=z(t_{f},e_{i})-z(t_{f},0)$ and
$J_{\varphi}(0):=\left[\beta_{1}(t)\ |\ \dots\ |\ \beta_{n}(t)\right]$.

Step 3

(Missing IC) Solve $J_{\varphi}(0)\,\lambda^{DR}_{0}:=-\varphi(0)$ for
$\lambda^{DR}_{0}$.

Step 4

(Projector onto ${\cal A}$) Solve (30) with ICs in (31)(iii) to find $x(t)$
and $\lambda^{DR}(t)$.
Set $P_{\cal{A}}(x^{-},u^{-})(t):=(y(t),v(t))$ where $y(t)=x(t)$ and
$v(t)=u^{-}(t)-B(t)^{T}\lambda^{DR}(t)$.

###  4.2 A conjecture for the costates for Problem (P)

Recall that the optimal control for Problem (P) is given by cases in (11). A
_junction time_ $\overline{t}_{j}$ is a time when the control $u_{j}(t)$ falls
into two cases of (11) simultaneously, i.e. a point in time where a control
constraint transitions from “active” to “inactive,” or vice versa. This
definition of a junction time becomes important in the following conjecture,
which has been formulated and tested by means of extensive numerical
experiments.

######  Conjecture 1.

Let $\lambda^{DR}(t)$ be the costate variable emerging from the projector into
$\cal{A}$ computed in Algorithm 2 and $\lambda(t)$ be the costate variable
emanating from Problem (P). Let $\overline{t}_{j}$ be a junction time for some
$u_{j}$, $j=1,\dots,m$, such that
$b_{j}(\overline{t}_{j})^{T}\lambda^{DR}(\overline{t}_{j})\neq 0$. Define

$\alpha:=\dfrac{-r_{j}(\overline{t}_{j})u_{j}(\overline{t}_{j})}{b_{j}(\overline{t}_{j})^{T}\lambda^{DR}(\overline{t}_{j})}\,.$

Then

$\lambda(t)=\alpha\lambda^{DR}(t)\,.$ (32)

######  Remark 4.

The ability to obtain the costate variable by Conjecture 1 is desirable as a
tool for checking that the necessary condition of optimality in (11) is
satisfied. Without this conjecture we are unable to verify whether the
optimality condition is satisfied when using the DR algorithm, except when a
dual version of the DR algorithm is employed, as in [13]. $\Box$

Once we have calculated $\lambda$ in this way we can also find a multiplier
$\mu^{k}$, $k=1,2$, numerically for the case when only one state constraint is
active at a given time. Suppose that only the $i$th state box constraint
becomes active. By rearranging Equation (1), using numerical differentiation
to find $\dot{\lambda}$ and assuming $\mu_{i}^{2}(t)=0$,

$\mu_{i}^{1}(t)=-Q(t)x(t)-A(t)^{T}\lambda(t)-\dot{\lambda}(t).$ (33)

If $\mu_{i}^{1}(t)=0$, then we compute

$\mu_{i}^{2}(t)=Q(t)x(t)+A(t)^{T}\lambda(t)+\dot{\lambda}(t).$ (34)

With (33), or with (34), the complementarity conditions in (3) or (4) can now
be checked numerically.

##  5 Numerical Experiments

We will now introduce two example problems. Along with posing the optimal
control problems we also give plots of their optimal controls, states,
costates and multipliers with vertical lines signifying the regions where the
state constraints become active.

###  5.1 Harmonic oscillator

Problem (PHO) below contains the dynamics of a harmonic oscillator which is
typically used to model a point mass with a spring. The dynamics are given as
$\ddot{y}(t)+\omega_{0}^{2}y(t)=f(t)$ where $\omega_{0}>0$ is known as the
natural frequency and $f(t)$ is some forcing. In a physical system $y$
represents the position of a unit mass, $\dot{y}$ is the velocity of said
mass, the natural frequency is expressed as $\omega_{0}=\sqrt{k}$ where $k$ is
the stiffness of the spring producing the harmonic motion and $f$ is the
restoring force. In addition to the restoring force we will introduce another
force $u_{1}$ that will affect the velocity $\dot{y}$ directly. We let
$x_{1}:=y$, $x_{2}:=\dot{y}$ and $u_{2}:=f$ to arrive at
$\dot{x}_{1}(t)=x_{2}(t)+u_{1}(t)$,
$\dot{x}_{2}(t)=-\omega_{0}^{2}x_{1}(t)+u_{2}(t)$.

In this example problem the objective contains the squared sum of all four
variables in the system. It is common to see this problem with the objective
of minimizing the energy of the control variable but in this case we have also
included the state variables to test the algorithm with a slightly more
involved objective. The focus of this research is control- and state-
constrained problems so the constraints are added as in Problem (P).

$\displaystyle\mbox{(PHO) }\left\\{\begin{array}[]{rl}\displaystyle\min&\ \
\displaystyle\frac{1}{2}\int_{0}^{2\pi}x_{1}^{2}(t)+x_{2}^{2}(t)+u_{1}^{2}(t)+u_{2}^{2}(t)\,dt\\\\[14.22636pt]
\mbox{subject to}&\ \ \dot{x}_{1}(t)=x_{2}(t)+u_{1}(t)\,,\ \ x_{1}(0)=0\,,\ \
x_{2}(2\pi)=0\,,\\\\[5.69054pt] &\ \ \dot{x}_{2}(t)=-4x_{1}(t)+u_{2}(t)\,,\ \
x_{2}(0)=1\,,\ \ x(2\pi)=0\,,\\\\[5.69054pt] &\ \ \underline{u}\leq
u(t)\leq\overline{u}\,,\ \ \underline{x}\leq
x(t)\leq\overline{x}\end{array}\right.$

###  5.2 Simple spring–mass system

The simple spring-mass system is another physical system that can be easily
visualised, see [30]. This problem contains two masses and two springs
connected in sequence with dynamics given by
$m_{1}\ddot{y}_{1}(t)+(k_{1}+k_{2})y_{1}(t)-k_{2}y_{2}=f_{1}(t)$,
$m_{2}\ddot{y}_{2}(t)-k_{2}y_{1}(t)+k_{2}y_{2}(t)=f_{2}(t)$ where
$m_{1},m_{2}$ are the two masses, $k_{1},k_{2}$ are the spring coefficients
(stiffness) and $f_{1}(t),f_{2}(t)$ are the forces applied to $m_{1},m_{2}$
respectively. Let $x_{1}:=y_{1}$, $x_{2}:=\dot{y}_{1}$, $x_{3}:=y_{2}$,
$x_{4}:=\dot{y}_{2}$, $u_{1}:=f_{1}$, $u_{2}:=f_{2}$, $m_{1}=m_{2}=1$,
$k_{1}=1$ and $k_{2}=2$ then we retrieve the system in Problem (PSM). This
dynamical system furnishes an optimal control problem with four state
variables and two control variables. As in (PHO) we add state and control
constraints and set the objective as the integral of the squared sum of the
state and control variables.

$\displaystyle\mbox{(PSM) }\left\\{\begin{array}[]{rl}\displaystyle\min&\ \
\displaystyle\frac{1}{2}\int_{0}^{2\pi}x_{1}^{2}(t)+x_{2}^{2}(t)+x_{3}^{2}(t)+x_{4}^{2}(t)+u_{1}^{2}(t)+u_{2}^{2}(t)\,dt\\\\[14.22636pt]
\mbox{subject to}&\ \ \dot{x}_{1}(t)=x_{2}(t)\,,\ \ x_{1}(0)=0\,,\ \
x_{1}(2\pi)=0\,,\\\\[5.69054pt] &\ \
\dot{x}_{2}(t)=-3x_{1}(t)+2x_{3}(t)+u_{1}(t)\,,\ \ x_{2}(0)=1\,,\ \
x_{2}(2\pi)=0\,,\\\\[5.69054pt] &\ \ \dot{x}_{3}(t)=x_{4}(t)\,,\ \
x_{3}(0)=1\,,\ \ x_{3}(2\pi)=0\,,\\\\[5.69054pt] &\ \
\dot{x}_{4}(t)=2x_{1}(t)-2x_{3}(t)+u_{2}(t)\,,\ \ x_{4}(0)=-1\,,\ \
x_{4}(2\pi)=0\,,\\\\[5.69054pt] &\ \ \underline{u}\leq
u(t)\leq\overline{u}\,,\ \ \underline{x}\leq
x(t)\leq\overline{x}\end{array}\right.$

###  5.3 Numerical discussion and comparisons

Technical Specifications. In the numerical experiments we used Matlab version
2023b with DR and compared with AMPL–Ipopt optimization computational suite
[19, 33] with Ipopt version 3.12.13. We chose to make comparisons to Ipopt
since it is a free and easily obtainable solver used for problems such as
those presented in this paper (also see [5]). All computations set
$\epsilon=10^{-8}$ from Algorithm 1 or in the case of Ipopt
$\text{tol}=10^{-8}$ and were run on a PC with an i5-10500T 2.30GHz processor
with 8GB RAM. For the two examples (PHO) and (PSM) we experimented with a case
that only had control constraints and a case that also had an added state
constraint. The problem specifications can be found in Table 1 along with the
choices of $\gamma$ that were used in the implementation of DR.

Problem | Case | $\gamma$ | Constraints
---|---|---|---
(PHO) | Case 1 | $0.60$ | $-0.4\leq u_{1}(t)\leq 0.1,\ -0.5\leq u_{2}(t)\leq 0.1$
Case 2 | $0.95$ | $-0.4\leq u_{1}(t)\leq 0.1,\ -0.5\leq u_{2}(t)\leq 0.1,\ -0.025\leq x_{1}(t)$
(PSM) | Case 1 | $0.55$ | $-0.5\leq u_{1}(t)\leq 0.5,\ -0.4\leq u_{2}(t)\leq 0.4$
Case 2 | $0.95$ | $-0.5\leq u_{1}(t)\leq 0.5,\ -0.4\leq u_{2}(t)\leq 0.4,\ -0.2\leq x_{1}(t)$
Table 1: Problem cases.

Generation of “True” Solutions. Since we do not have analytical solutions to
these problems we have generated higher accuracy numerical solutions to our
problems that we will use as “true” solutions in our error calculations. In
the Case 1 examples (only control constraints present) DR was able to
successfully converge to an acceptable solution using
$N=10^{7},\epsilon=10^{-12}$ without reaching 200 iterations (the maximum
number of iterations we allowed). The 8GB RAM provided insufficient memory for
Ipopt to generate a solution with such a high number of discretization points
so DR was used to generate the “true” solutions in this case. In Case 2 (state
constrained case) DR was unable to converge to an accepted solution in less
than 200 iterations so we instead relied on Ipopt. Due to the memory
limitations the “true” solution for (PHO) Case 2 was generated with
$N=10^{6},\ \epsilon=10^{-12}$ and (PSM) Case 2 used $N=7\times 10^{5},\
\epsilon=10^{-12}$.

Choice of an Algorithmic Parameter. The values of $\gamma$ in Table 1 were
decided by generating plots displaying the number of iterations required for
DR to find an acceptable solution for 500 values of $\gamma$ in the interval
$(0,1)$. For both (PHO) and (PSM) in Case 2 for all values of $\gamma$ that
were tested DR required more than the maximum 200 iterations to converge. For
these cases we instead generated plots that compared the errors in the states
and controls for the different values of $\gamma$. A specific value of
$\gamma$ that would provide the smallest errors was not obvious since many
values provided similar performance but values closer to 1 appeared optimal
thus the choice of $\gamma=0.95$ was made for these experiments.

Figure 1: (PHO) Case 2 plots (see Table 1) using DR with $N=10^{3},-0.025\leq
x_{1}(t)$. Vertical lines indicate the interval in which the state constraint
becomes active.

Figure 2: (PSM) Case 2 plots (see Table 1) using DR with $N=10^{3},-0.2\leq
x_{1}(t)$. Vertical lines indicate the interval in which the state constraint
becomes active. Figure 3: Multipliers $\mu_{1}^{2},\mu_{2}^{2}$ for Case 2
(PHO) with $N=10^{3}$ using DR (solid lines) and Ipopt (dotted lines). Note
that $\mu_{1}^{2}(t)$ found by Ipopt attains a maximum value of 33, which is
not shown in the graph.

Graphical Results. In Figures 1–2 we have plots for (PHO) and (PSM) using DR.
We have generated similar figures using Ipopt as well but since they mostly
overlap Figures 1–2 we will not show them all here but will point out the
differences. In Figure 3 we give the multiplier vector $\mu^{2}$ components
for Case 2 (PHO) using DR and Ipopt. As indicated by the black vertical lines
in the bottom left plot of Figure 1 when using DR the constraint on $x_{1}$ is
active over a time interval that aligns with the interval where $\mu_{1}^{2}$
is positive, as expected from Equations (3)–(4).

In the results from Ipopt we observe that the state constraint only became
active at a single point in time but Figure 3 shows that $\mu_{1}^{2}$ from
Ipopt (yellow dotted line) is positive for a larger interval on time thus
Equation (4) is violated. This appears to be a numerical error that is not
present in the DR results since there is an interval of points around the
spike where the state constraint $-0.025\leq x_{1}(t)$ becomes active. When
$N=10^{4},10^{5}$ for (PHO) we observe that when using Ipopt the lower bound
on $x_{1}$ is never reached though again we see the interval of points that
are almost equal to the lower bound.

In Figure 3 we also see a discrepancy in $\mu_{2}^{2}$. Since we have not
imposed a constraint on the variable $x_{2}$ Equation (4) implies that
$\mu_{2}^{2}(t)=0$ for all $t\in[t_{0},t_{f}]$. We can see in Figure 3 that
$\mu_{2}^{2}$ from Ipopt (purple dotted line) fails to satisfy this
requirement while $\mu_{2}^{2}$ from DR is, at least to the eye, equal to
zero.

Another difference between the multipliers $\mu$ from DR and Ipopt is the
maximum values reached by the functions. In Figure 3 we see that $\mu_{1}^{2}$
obtained via DR and Ipopt have similar shapes in their plots but the maximum
value reached using DR is approximately $0.7$ while from Ipopt the maximum
value is approximately $33$. For (PSM) Case 2 with $N=10^{3}$ the maximum
value that $\mu_{1}^{2}$ obtained using DR was approximately 16 while Ipopt
reached approximately 310. Along with the functions having very different
maximum values we noted that when generating these plots for
$N=10^{3},10^{4},10^{5}$ the results from DR were clearly converging to a
single function $\mu_{1}^{2}$ while this was not obvious from the Ipopt
results. For (PSM) we see the approximate maximum values of $\mu_{1}^{2}$
obtained when using Ipopt were 310 for $N=10^{3}$, 2500 for $N=10^{4}$ and
3000 for $N=10^{5}$. In the same example using DR the maximum values of the
function were 16 for $N=10^{3}$, 14 for $N=10^{4}$ and 14 for $N=10^{5}$.

We also see some slight variation between the plots from DR and Ipopt at the
junction times in the control variables (not shown in this paper in order to
avoid excessive amount of visual material). The intervals where the control
variables attain their boundaries using DR always appear slightly larger than
those from Ipopt. The control variables in the region where they transition
between active and inactive constraints appear more rounded at these corners
when using Ipopt and have a more sharp transition when using DR.

Errors and CPU Times. Table 2 contains the errors in the controls, states,
costates for DR and Ipopt while Table 3 contains the errors in the
multipliers, objective values and CPU times. The CPU times were computed as
averages of over 200 runs (up to 1,000 runs in the faster examples). The
values within boxes are the smaller errors and CPU times between DR and Ipopt.
At a glance we can see that more often than not DR produced smaller errors and
faster CPU times when compared with Ipopt. Upon closer inspection we see that
in many of the Case 2 results the errors from DR and Ipopt are comparable. We
note that the ”true” solutions we are using to calculate these errors are
those explained earlier in this subsection except for the multipliers
$\mu^{2}$ in the case 1 examples. In those examples we take the zero vector as
our “true” solution.

We observe that most of the results in Table 2 show a smaller error in the
control variable from DR, especially in the case 1 examples where there are no
state constraints. In the state variables we see that Ipopt has the smaller
errors when $N=10^{3},10^{4}$ though there is little difference compared with
DR and we see an improvement in DR when $N=10^{5}$. Like with the error in the
control variables we see that the error in the costates is smaller for DR in
almost all examples.

From Table 3 the errors in the multipliers show an improvement in DR compared
with Ipopt though the “true” solution in the case 2 examples was obtained
using results from Ipopt which as previously mentioned did not appear to
converge to a specific value as we increased $N$ so the quality of this “true”
solution is not guaranteed. DR produced slightly smaller objective values that
were closest to the “true” solution in almost all cases though the difference
compared with Ipopt is marginal. We see that the CPU times are faster for DR
especially in the examples where $N=10^{5}$.

|  |  $L^{\infty}$ error in control |  $L^{\infty}$ error in states |  $L^{\infty}$ error in costates
---|---|---|---|---
$N$ | Problem | DR | Ipopt | DR | Ipopt | DR | Ipopt
$10^{3}$ | (PHO) - Case 1 | $\boxed{7.9\times 10^{-3}}$ | $1.4\times 10^{-2}$ | $6.9\times 10^{-3}$ | $\boxed{2.7\times 10^{-3}}$ | $\boxed{6.3\times 10^{-3}}$ | $1.4\times 10^{-2}$
| (PHO) - Case 2 | $\boxed{1.4\times 10^{-2}}$ | $1.6\times 10^{-2}$ | $4.5\times 10^{-3}$ | $\boxed{2.9\times 10^{-3}}$ | $\boxed{1.4\times 10^{-2}}$ | $1.6\times 10^{-2}$
| (PSM) - Case 1 | $\boxed{2.3\times 10^{-2}}$ | $3.8\times 10^{-2}$ | $2.0\times 10^{-2}$ | $\boxed{1.8\times 10^{-2}}$ | $\boxed{5.5\times 10^{-2}}$ | $1.5\times 10^{-1}$
| (PSM) - Case 2 | $\boxed{7.1\times 10^{-2}}$ | $1.1\times 10^{-1}$ | $3.8\times 10^{-1}$ | $\boxed{3.7\times 10^{-1}}$ | $\boxed{8.4\times 10^{-1}}$ | $1.5\times 10^{0}$
$10^{4}$ | (PHO) - Case 1 | $\boxed{7.8\times 10^{-4}}$ | $3.4\times 10^{-3}$ | $6.7\times 10^{-4}$ | $\boxed{3.6\times 10^{-4}}$ | $\boxed{7.5\times 10^{-4}}$ | $1.6\times 10^{-3}$
| (PHO) - Case 2 | $\boxed{1.3\times 10^{-3}}$ | $3.0\times 10^{-3}$ | $4.7\times 10^{-4}$ | $\boxed{4.2\times 10^{-4}}$ | $5.0\times 10^{-3}$ | $\boxed{4.8\times 10^{-3}}$
| (PSM) - Case 1 | $\boxed{2.2\times 10^{-3}}$ | $3.6\times 10^{-3}$ | $2.0\times 10^{-3}$ | $\boxed{1.8\times 10^{-3}}$ | $\boxed{4.8\times 10^{-3}}$ | $1.5\times 10^{-2}$
| (PSM) - Case 2 | $2.4\times 10^{-2}$ | $\boxed{1.1\times 10^{-2}}$ | $\boxed{3.7\times 10^{-1}}$ | $\boxed{3.7\times 10^{-1}}$ | $\boxed{7.9\times 10^{-1}}$ | $1.5\times 10^{0}$
$10^{5}$ | (PHO) - Case 1 | $\boxed{7.7\times 10^{-5}}$ | $6.3\times 10^{-3}$ | $\boxed{6.7\times 10^{-5}}$ | $9.5\times 10^{-4}$ | $\boxed{6.5\times 10^{-5}}$ | $2.1\times 10^{-3}$
| (PHO) - Case 2 | $\boxed{7.8\times 10^{-4}}$ | $7.5\times 10^{-3}$ | $\boxed{6.4\times 10^{-5}}$ | $1.4\times 10^{-3}$ | $\boxed{5.7\times 10^{-3}}$ | $7.3\times 10^{-3}$
| (PSM) - Case 1 | $\boxed{2.2\times 10^{-4}}$ | $1.1\times 10^{-2}$ | $\boxed{2.0\times 10^{-4}}$ | $6.8\times 10^{-4}$ | $\boxed{4.3\times 10^{-4}}$ | $4.0\times 10^{-3}$
| (PSM) - Case 2 | $2.6\times 10^{-2}$ | $\boxed{6.0\times 10^{-3}}$ | $\boxed{3.7\times 10^{-1}}$ | $\boxed{3.7\times 10^{-1}}$ | $8.0\times 10^{-1}$ | $\boxed{7.9\times 10^{-1}}$
Table 2: Errors in controls $u$, states $x$ and costates $\lambda$ and CPU
times for the DR algorithm and AMPL–Ipopt,
with $\varepsilon=10^{-8}$ and specifications from Table 1. |  |  $L^{\infty}$ error in $\mu$ |  $L^{\infty}$ error in objective values | CPU time [sec]
---|---|---|---|---
$N$ | Problem | DR | Ipopt | DR | Ipopt | DR | Ipopt
$10^{3}$ | (PHO) - Case 1 | $\boxed{4.9\times 10^{-3}}$ | $3.1\times 10^{-2}$ | $\boxed{2.9\times 10^{-3}}$ | $4.8\times 10^{-3}$ | $\boxed{4.8\times 10^{-2}}$ | $2.9\times 10^{-1}$
| (PHO) - Case 2 | $\boxed{5.4\times 10^{-1}}$ | $1.5\times 10^{0}$ | $\boxed{2.9\times 10^{-3}}$ | $4.9\times 10^{-3}$ | $\boxed{3.3\times 10^{-1}}$ | $3.7\times 10^{-1}$
| (PSM) - Case 1 | $\boxed{2.4\times 10^{-3}}$ | $8.6\times 10^{-2}$ | $\boxed{4.8\times 10^{-2}}$ | $5.1\times 10^{-2}$ | $\boxed{1.0\times 10^{-1}}$ | $4.4\times 10^{-1}$
| (PSM) - Case 2 | $1.6\times 10^{1}$ | $\boxed{3.1\times 10^{0}}$ | $\boxed{6.8\times 10^{-2}}$ | $7.4\times 10^{-2}$ | $\boxed{5.2\times 10^{-1}}$ | $5.3\times 10^{-1}$
$10^{4}$ | (PHO) - Case 1 | $\boxed{1.5\times 10^{-4}}$ | $3.0\times 10^{-3}$ | $\boxed{2.8\times 10^{-4}}$ | $4.8\times 10^{-4}$ | $\boxed{4.7\times 10^{-1}}$ | $2.7\times 10^{0}$
| (PHO) - Case 2 | $\boxed{1.7\times 10^{1}}$ | $\boxed{1.7\times 10^{1}}$ | $\boxed{2.8\times 10^{-4}}$ | $4.9\times 10^{-4}$ | $3.5\times 10^{0}$ | $\boxed{3.3\times 10^{0}}$
| (PSM) - Case 1 | $\boxed{7.1\times 10^{-4}}$ | $8.2\times 10^{-3}$ | $\boxed{4.6\times 10^{-3}}$ | $5.0\times 10^{-3}$ | $\boxed{1.1\times 10^{0}}$ | $4.6\times 10^{0}$
| (PSM) - Case 2 | $1.4\times 10^{1}$ | $\boxed{2.5\times 10^{0}}$ | $\boxed{4.4\times 10^{-3}}$ | $6.9\times 10^{-3}$ | $\boxed{5.9\times 10^{0}}$ | $6.4\times 10^{0}$
$10^{5}$ | (PHO) - Case 1 | $\boxed{3.7\times 10^{-5}}$ | $3.0\times 10^{-4}$ | $\boxed{2.8\times 10^{-5}}$ | $7.7\times 10^{-5}$ | $\boxed{5.1\times 10^{0}}$ | $2.6\times 10^{1}$
| (PHO) - Case 2 | $\boxed{1.8\times 10^{1}}$ | $\boxed{1.8\times 10^{1}}$ | $\boxed{2.4\times 10^{-5}}$ | $1.4\times 10^{-4}$ | $3.7\times 10^{1}$ | $\boxed{3.1\times 10^{1}}$
| (PSM) - Case 1 | $\boxed{1.1\times 10^{-4}}$ | $8.2\times 10^{-4}$ | $\boxed{4.5\times 10^{-4}}$ | $5.3\times 10^{-4}$ | $\boxed{1.2\times 10^{1}}$ | $1.2\times 10^{2}$
| (PSM) - Case 2 | $\boxed{2.2\times 10^{2}}$ | $3.7\times 10^{3}$ | $1.5\times 10^{-3}$ | $\boxed{7.2\times 10^{-4}}$ | $\boxed{6.6\times 10^{1}}$ | $1.6\times 10^{2}$
Table 3: Errors in multipliers $\mu$, errors in objective values and CPU times
for the DR algorithm and AMPL–Ipopt,
with $\varepsilon=10^{-8}$ and specifications from Table 1.

##  6 Conclusion

We have applied the Douglas–Rachford (DR) algorithm to find a numerical
solution of LQ control problems with state and control constraints, after re-
formulating these problems as the sum of two convex functions and deriving
expressions for the proximal mappings of these functions (Theorems 1–2). These
proximal mappings were used in the DR iterations. Within the DR algorithm
(Algorithm 1), we proposed a procedure (Algorithm 2) for finding the
projection onto the affine set defined by the ODEs numerically.

We carried out extensive numerical experiments on two nontrivial example
problems and illustrated both the working of the algorithm and its efficiency
(in both accuracy and speed) compared with the traditional approach of direct
discretization. We observed that in general the DR algorithm produced smaller
errors and faster run times for these problems most notably in the examples
where we have increased the number of discretization points. From these
numerical results the DR algorithm could in general be recommended over Ipopt
when generating high quality solutions.

Based on further extensive experiments, we conjectured on how the costate
variables can be determined. We successfully used the costate variables
constructed as in the conjecture, as well as the state constraint multipliers
calculated using these costate variables, for the numerical verification of
the optimality conditions.

We recall that Algorithm 2 involves repeated numerical solution of the ODEs in
(30) with various initial conditions. For solving (30) we implemented
(explicit) Euler’s method which requires only a continuous right-hand side of
the ODEs. Algorithm 2 appears to be successful for the worked examples partly
owing to the fact that in these examples the optimal control is continuous. We
tried to use our approach for the machine tool manipulator problem from [16]
which has 7 state variables, one control variable, and upper and lower bounds
imposed on the single control variable and one of the state variables.
However, our approach did not seem to yield a solution (so far) for this
problem, conceivably owing to the fact that the optimal control variable, as
well as the optimal costate variable vector, is not continuous, in that these
variables have jumps a number of times during the process—see [16, Figure 5].
Note that discontinuities in the control make the right-hand side of (30)
discontinuous rendering Euler’s method ineffective. Therefore, problems of the
kind in [16] require further investigation.

We believe that our approach can be extended to more general convex optimal
control problems, for example those with a non-quadratic objective function or
mixed state and control constraints, as long as the pertaining proximal
operators are not too challenging to derive.

It would also be interesting to employ and test, in the future, other
projection type methods, for example the Peaceman–Rachford algorithm [8,
Section 26.4 and Proposition 28.8], which, to the knowledge of the authors,
has not been applied to optimal control problems.

## Acknowledgments

BIC was supported by an Australian Government Research Training Program
Scholarship. No funding was received by RSB and CYK to assist with the
preparation of this manuscript.

## Conflict of Interest

The authors have no competing, or conflict of, interests to declare that are
relevant to the content of this article.

## Use of AI Tools Declaration

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

## References

  * [1] H. M. Amman, D. A. Kendrick, Computing the steady state of linear quadratic optimization models with rational expectations. Econ. Lett., 58(2), 185–191, 1998.
  * [2] F. J. Aragón Artacho, J. M. Borwein and M. K. Tam, Douglas–Rachford feasibility methods For matrix completion problems. ANZIAM Journal, 55(4), 299–326, 2014.
  * [3] F. J. Aragón Artacho, R. Campoy and V. Elser, An enhanced formulation for solving graph coloring problems with the Douglas–Rachford algorithm. J. Glob. Optim., 77(2), 383–403, 2020.
  * [4] U. M. Ascher, R. M. M. Mattheij, and R. D. Russell, Numerical Solution of Boundary Value Problems for Ordinary Differential Equations. SIAM Publications, Philadelphia, 1995.
  * [5] N. Banihashemi and C. Y. Kaya, Inexact restoration for Euler discretization of box-constrained optimal control problems. J. Optim. Theory Appl., 156 (2013), 726–760.
  * [6] H. H. Bauschke, 8 Queens, Sudoku, and Projection Methods. https://carma.newcastle.edu.au/resources/jon/Preprints/Books/CUP/Material/Lions-Mercier/Heinz_Bauschke.pdf, 2008.
  * [7] H. H. Bauschke, R. S. Burachik, and C. Y. Kaya, Constraint splitting and projection methods for optimal control of double integrator in Splitting Algorithms, Modern Operator Theory, and Applications. Springer, 45–68, 2019.
  * [8] H. H. Bauschke and P. L. Combettes, Convex Analysis and Monotone Operator Theory in Hilbert Spaces, Second edition. Springer, 2017.
  * [9] H. H. Bauschke and V. R. Koch, Projection methods: Swiss Army knives for solving feasibility and best approximation problems with halfspaces. Contemporary Mathematics, 636, 1–40, 2015.
  * [10] H. H. Bauschke and W. M. Moursi, On the Douglas–Rachford algorithm. Math. Program., Ser. A, 164, 263–284, 2017.
  * [11] R. S. Burachik, B. I. Caldwell, C. Y. Kaya, Projection methods for control-constrained minimum-energy control problems. arXiv:2210.17279v1. https://arxiv.org/abs/2210.17279, 2022.
  * [12] R. S. Burachik, B. I. Caldwell, C. Y. Kaya, Douglas–Rachford algorithm for control-constrained minimum-energy control problems. To appear in ESAIM Control Optim. Calc. Var. arXiv:2210.17279v2, https://arxiv.org/abs/….., 2024.
  * [13] R. S. Burachik, B. I. Caldwell, C. Y. Kaya, W. M. Moursi, Optimal control duality and the Douglas–Rachford algorithm. To appear in SIAM J. Control Optim.
https://doi.org/10.48550/arXiv.2303.06527, 2023.

  * [14] R. S. Burachik, C. Y. Kaya, and S. N. Majeed, A duality approach for solving control-constrained linear-quadratic optimal control problems. SIAM J. Control Optim., 52 (2014), 1771–1782.
  * [15] C. Büskens, H. Maurer, SQP-methods for solving optimal control problems with control and state constraints: Adjoint variables, sensitivity analysis and real-time control. J. Comput. Appl. Math., 120(1), 85–108, 2000.
  * [16] B. Christiansen, H. Maurer, and O. Zirn, Optimal control of machine tool manipulators. Diehl M, Glineur F, Jarlebring E, Michiels W (eds), Recent Advances in Optimization and Its Applications in Engineering. Springer, Berlin, Heidelberg, 451–460, 2010.
  * [17] J. Douglas and H. H. Rachford, On the numerical solution of heat conduction problems in two and three space variables. Trans. Amer. Math. Soc., 82, 421–439, 1956.
  * [18] J. Eckstein and D. P. Bertsekas, On the Douglas-Rachford splitting method and the proximal point algorithm for maximal monotone operators. Math. Progr., Ser. A, 55, 293–318, 1992.
  * [19] R. Fourer, D. M. Gay, and B. W. Kernighan, AMPL: A Modeling Language for Mathematical Programming, 2nd edition. Brooks/Cole Publishing Company / Cengage Learning, 2003.
  * [20] S. Gravel and V. Elser, Divide and concur: A general approach to constraint satisfaction. Physical Review E, Statistical, Nonlinear, and Soft Matter Physics, 78(3), 036706–036706, 2008.
  * [21] R. F. Hartl, S. P. Sethi, R. G. Vickson, A survey of the maximum principles for optimal control problems with state constraints. SIAM Review, 52(2), 181–218, 1995.
  * [22] C. Y. Kaya, Inexact Restoration for Runge-Kutta discretization of optimal control problems, SIAM J. Numer. Anal., 48(4), 1492–1517, 2010.
  * [23] H. B. Keller, Numerical Methods for Two-Point Boundary-Value Problems. Blaisdell, London, 1968.
  * [24] B. Kugelmann, H. J. Pesch, New general guidance method in constrained optimal control. I, Numerical method. J. Optim. Theory Appl., 67(3), 421–435, 1990.
  * [25] P.-L. Lions and B. Mercier, Splitting algorithms for the sum of two nonlinear operators. SIAM J. Numer. Anal., 16, 964–979, 1979.
  * [26] H. Maurer, J.-H.R. Kim and G. Vossen, On a state-constrained control problem in optimal production and maintenance. Optimal Control and Dynamic Games, Applications in Finance, Management Science and Economics, (C. Deissenberg, R.F. Hartl, eds.), Springer Verlag, 289–308, 2005.
  * [27] H. Maurer, H. J. Oberle, Second order sufficient conditions for optimal control problems with free final time: The Riccati approach. SIAM J. Control Optim., 41 (2), 380–403, 2003.
  * [28] T. Mouktonglang, Innate immune response via perturbed LQ-control problem. Advanced Studies in Biology, 3, 327–332, 2011.
  * [29] W. J. Rugh, Linear System Theory, 2nd edition, Pearson, 1996.
  * [30] T. L. Schmitz and K. S. Smith, Mechanical Vibrations Modeling and Measurement, 2nd edition. Springer, 2021.
  * [31] J. Stoer and R. Bulirsch, Introduction to Numerical Analysis, 2nd edition. Springer-Verlag, New York, 1993.
  * [32] B. F. Svaiter, On weak convergence of the Douglas-Rachford method. SIAM J. Control Optim., 49, 280–287, 2011.
  * [33] A. Wächter and L. T. Biegler, On the implementation of a primal-dual interior point filter line search algorithm for large-scale nonlinear programming. Math. Progr., 106, 25–57, 2006.


# Approximation by Nörlund means with respect to Vilenkin system in Lebesgue
spaces

N. Anakidze, N. Areshidze, L. Baramidze  N.Anakidze, The University of
Georgia, School of Science and Technology, 77a Merab Kostava St, Tbilisi,
0171, Georgia<EMAIL_ADDRESS>N.Areshidze, Tbilisi State University,
Faculty of Exact and Natural Sciences, Department of Mathematics, Chavchavadze
str. 1, tbilisi 0128, georgia<EMAIL_ADDRESS>L.Baramidze, Tbilisi
State University, Faculty of Exact and Natural Sciences, Department of
Mathematics, Chavchavadze str. 1, tbilisi 0128, georgia<EMAIL_ADDRESS>
###### Abstract.

In this paper we improve and complement a result by Móricz and Siddiqi [19].
In particular, we prove that their estimate of the Nörlund means with respect
to the Vilenkin system holds also without their additional condition.
Moreover, we prove a similar approximation result in Lebesgue spaces for any
$1\leq p<\infty$.

2010 Mathematics Subject Classification: 42C10, 42B30.

Key words and phrases: Vilenkin group, Vilenkin system, Fejér means, Nörlund
means, approximation.

##  1\. Introduction

Concerning some definitions and notations used in this introduction we refer
to Section 2.

It is well-known (see e.g. [16], [28] and [39]) that, for any $1\leq
p\leq\infty$ and $f\in L_{p}(G_{m}),$ there exists an absolute constant
$C_{p},$ depending only on $p$ such that

$\left\|\sigma_{n}f\right\|_{p}\leq C_{p}\left\|f\right\|_{p}.$

Moreover, (for details see [28]) if $1\leq p\leq\infty$, $M_{N}\leq
n<M_{N+1}$, $f\in L^{p}(G_{m})$ and $n\in\mathbb{N},$ then

(1) $\displaystyle\left\|\sigma_{n}f-f\right\|_{p}\leq
R^{2}\sum_{s=0}^{N}\frac{M_{s}}{M_{N}}\omega_{p}\left(1/M_{s},f\right),$

where $R:=\sup_{k\in\mathbb{N}}m_{k}$ and $\omega_{p}(\delta,f)$ is the
modulus of continuity of $L^{p},1\leq p\leq\infty$ functions defined by

$\displaystyle\omega_{p}(\delta,f)=\sup_{|t|<\delta}\|f(x+t)-f(x)\|,\quad\quad\delta>0.$

It follows that if $f\in lip\left(\alpha,p\right),$ i. e.,

$\displaystyle Lip(\alpha,p):=\\{f\in
L^{p}:\omega_{p}(\delta,f)=O(\delta^{\alpha})\quad\textrm{as}\quad\delta\to
0\\},$

then

$\left\|\sigma_{n}f-f\right\|_{p}=\left\\{\begin{array}[]{c}O\left(1/M_{N}\right),\text{
\ \ \ if \ \ \ }\alpha>1,\\\ O\left(N/M_{N}\right),\text{ \ \ \ if \ \ \
}\alpha=1,\\\ O\left(1/M_{n}^{\alpha}\right),\text{ \ \ \ \ \ \ \ if
}\alpha<1.\end{array}\right.$

Moreover, (for details see [28]) if $1\leq p<\infty,$ $f\in L^{p}(G)$ and

$\left\|\sigma_{M_{n}}f-f\right\|_{p}=o\left(1/M_{n}\right),\ \text{as}\
n\rightarrow\infty,$

then $f$ is a constant function.

The weak-$(1,1)$ type inequality for the maximal operators of Vilenkin-Fejer
means $\sigma^{*}$, defined by

$\sigma^{*}f=\sup_{n\in\mathbb{N}}|\sigma_{n}f|$

can be found in Schipp [29] for Walsh series and in Pál, Simon [27] and Weisz
[37] for bounded Vilenkin series. Boundedness of the maximal operators of
Vilenkin-Féjer means of the one- and two-dimensional cases can be found in
Fridli [11], Gát [13], Goginava [15], Nagy and Tephnadze [23, 24, 25, 26],
Simon [31, 32], Tutberidze [33], Weisz [38].

Convergence and summability of Nörlund means with respect to Vilenkin systems
were studied by Areshidze and Tephnadze [2], Blahota, Persson and Tephnadze
[8] (see also [3, 4, 6, 7]), Fridli, Manchanda and Siddiqi [12], Goginava
[14], Nagy [20, 21, 22] (see also [9] and [10]), Memic [17].

Móricz and Siddiqi [19] investigated the approximation properties of some
special Nörlund means of Walsh-Fourier series of $L^{p}$ functions in norm. In
particular, they proved that if $f\in L^{p}(G),$ $1\leq p\leq\infty,$
$n=2^{j}+k,$ $1\leq k\leq 2^{j}\ (n\in\mathbb{N}_{+})$ and
$(q_{k},k\in\mathbb{N})$ is a sequence of non-negative numbers, such that

(2) $\frac{n^{\gamma-1}}{Q_{n}^{\gamma}}\sum_{k=0}^{n-1}q^{\gamma}_{k}=O(1),\
\ \text{for some}\ \ 1<\gamma\leq 2,$

then, there exists an absolute constant $C_{p},$ depending only on $p$ such
that

(3)
$\|t_{n}f-f\|_{p}\leq\frac{C_{p}}{Q_{n}}\sum_{i=0}^{j-1}2^{i}q_{n-2^{i}}\omega_{p}\left(\frac{1}{2^{i}},f\right)+C_{p}\omega_{p}\left(\frac{1}{2^{j}},f\right),$

when the sequence $(q_{k},k\in\mathbb{N})$ is non-decreasing, while

$\|t_{n}f-f\|_{p}\leq\frac{C_{p}}{Q_{n}}\sum_{i=0}^{j-1}\left(Q_{n-2^{i}+1}-Q_{n-2^{i+1}+1}\right)\omega_{p}\left(\frac{1}{2^{i}},f\right)+C_{p}\omega_{p}\left(\frac{1}{2^{j}},f\right),$

when the sequence $(q_{k},k\in\mathbb{N})$ is non-increasing.

In this paper we improve and complement a result by Móricz and Siddiqi [19].
In particular, we prove that their estimate of the Nörlund means with respect
to the Vilenkin system holds also without their additional condition.
Moreover, we prove a similar approximation result in Lebesgue spaces for any
$1\leq p<\infty$.

##  2\. Preliminaries

Let $\mathbb{N}_{+}$ denote the set of the positive integers,
$\mathbb{N}:=\mathbb{N}_{+}\cup\\{0\\}.$ Let $m=:(m_{0},m_{1},...)$ be a
sequence of positive integers not less than 2. Denote by

$Z_{m_{k}}:=\\{0,1,...,m_{k}-1\\}$

the additive group of integers modulo $m_{k}.$ Define the group $G_{m}$ as the
complete direct product of the group $Z_{m_{k}}$ with the product of the
discrete topologies of $Z_{m_{k}}\textrm{'s}$.

The direct product $\mu$ of the measures

$\mu_{k}\left(\\{j\\}\right):=1/m_{k}\text{ \ }(j\in Z_{m_{k}})$

is the Haar measure on $G_{m}$ with $\mu\left(G_{m}\right)=1.$

If $\sup_{k\in\mathbb{N}}m_{k}<+\infty$, then we call $G_{m}$ a bounded
Vilenkin group. If $\\{m_{k}\\}_{k\geq 0}$ sequance is unbounded, then $G_{m}$
is said to be unbounded Vilenkin group. In this paper we consider only bounded
Vilenkin groups.

The elements of $G_{m}$ are represented by the sequences

$x:=(x_{0},x_{1},\dots,x_{k},\dots)\qquad\left(\text{ }x_{k}\in
Z_{m_{k}}\right).$

It is easy to give a base for the neighborhood of $G_{m}$, namely

$\displaystyle I_{0}\left(x\right):=G_{m},\ \ \ I_{n}(x):=\\{y\in G_{m}\mid
y_{0}=x_{0},\dots,y_{n-1}=x_{n-1}\\}\text{ }(x\in G_{m},\text{
}n\in\mathbb{N}).$

For the simplicity we also define by $I_{n}$ as $I_{n}:=I_{n}(0).$

Let us define a generalized number system based on $m$ in the following way:

$M_{0}=:1,\quad M_{k+1}=:m_{k}M_{k}\quad(k\in\mathbb{N})$

Then every $n\in\mathbb{N}$ can be uniquely expressed as

$n=\sum_{k=0}^{\infty}n_{j}M_{j},\ \ \ \text{where }\ \ \ n_{j}\in Z_{m_{j}}\
\ \ (j\in\mathbb{N})$

and only a finite number of $n_{j}`$s differ from zero. Let

$|n|=:\max\\{j\in\mathbb{N},n_{j}\neq 0\\}.$

In 1947 Vilenkin [34, 35, 36] investigated a group $G_{m}$ and introduced the
Vilenkin systems $\\{{\psi}_{j}\\}_{j=0}^{\infty}$ as

$\psi_{n}\left(x\right):=\prod_{k=0}^{\infty}r_{k}^{n_{k}}\left(x\right)\text{
\quad}\left(n\in\mathbb{N}\right).$

where $r_{k}(x)$ are the generalized Rademacher functions defined by

$r_{k}(x):=\exp(2\pi ix_{k}/m_{k}),\text{ \ \ }\left(k\in\mathbb{N}\right).$

These systems include as a special case the Walsh system when $m_{k}=2$ for
any $k\in\mathbb{N}$.

The norms (or quasi-norms) of Lebesgue spaces $L_{p}(G_{m})$ are defined by

$\left\|f\right\|_{p}^{p}:=\int_{G_{m}}\left|f\right|^{p}d\mu.$

The Vilenkin system is orthonormal and complete in $L^{2}\left(G_{m}\right)\,$
(for details see e.g. [1] and [30]).

If $f\in L^{1}\left(G_{m}\right)$, we can define the Fourier coefficients, the
partial sums of the Fourier series, the Fejér means, the Dirichlet and Fejér
kernels with respect to the Vilenkin system in the usual manner:

$\displaystyle\widehat{f}\left(k\right)$ $\displaystyle:$
$\displaystyle=\int_{G_{m}}f\overline{\psi}_{k}d\mu,\,\text{\quad}\left(k\in\mathbb{N}\right),$
$\displaystyle S_{n}f$ $\displaystyle:$
$\displaystyle=\sum_{k=0}^{n-1}\widehat{f}\left(k\right)\psi_{k},\text{
\quad}\left(n\in\mathbb{N}_{+},\text{ }S_{0}f:=0\right),$
$\displaystyle\sigma_{n}f$ $\displaystyle:$
$\displaystyle=\frac{1}{n}\sum_{k=1}^{n}S_{k}f,\text{ \quad\ \
}\left(n\in\mathbb{N}_{+}\right).$ $\displaystyle D_{n}$ $\displaystyle:$
$\displaystyle=\sum_{k=0}^{n-1}\psi_{k},\text{ \quad\ \
}\left(n\in\mathbb{N}_{+}\right).$ $\displaystyle K_{n}$ $\displaystyle:$
$\displaystyle=\frac{1}{n}\sum_{k=1}^{n}D_{k},\text{ \quad\ \
}\left(n\in\mathbb{N}_{+}\right).\ $

Recall that (for details see e.g. [1] and [28]),

(4) $\quad D_{M_{n}}\left(x\right)=\left\\{\begin{array}[]{l}\text{
}M_{n},\text{\thinspace\thinspace\thinspace if\thinspace\thinspace}x\in
I_{n},\\\ \text{
}0,\text{\thinspace\thinspace\thinspace\thinspace\thinspace\thinspace\thinspace\thinspace
if \thinspace\thinspace}x\notin I_{n},\end{array}\right.$

$\displaystyle D_{M_{n}-j}\left(x\right)$ $\displaystyle=$ $\displaystyle
D_{M_{n}}\left(x\right)-\overline{\psi}_{M_{n}-1}(-x)D_{j}(-x)$
$\displaystyle=$ $\displaystyle
D_{M_{n}}\left(x\right)-\psi_{M_{n}-1}(x)\overline{D_{j}}(x),\,\,\ 0\leq
j<M_{n}.$ (6) $n\left|K_{n}\right|\leq
2\sum_{l=0}^{|n|}M_{l}\left|K_{M_{l}}\right|,$

and

(7) $\displaystyle\int_{G_{m}}K_{n}(x)d\mu(x)=1,\ \ \ \ \
\sup_{n\in\mathbb{N}}\int_{G_{m}}\left|K_{n}(x)\right|d\mu(x)\leq 2.$

Moreover, if $n>t,$ $t,n\in\mathbb{N},$ then

(8)
$K_{M_{n}}\left(x\right)=\left\\{\begin{array}[]{ll}\frac{M_{t}}{1-r_{t}(x)},&x\in
I_{t}\backslash I_{t+1},\quad x-x_{t}e_{t}\in I_{n},\\\
\frac{M_{n}+1}{2},&x\in I_{n},\\\ 0,&\text{otherwise.}\end{array}\right.$

The $n$-th Nörlund mean $t_{n}$ of the Vilenkin-Fourier series of a integrable
function $f$ is defined by

(9) $t_{n}f:=\frac{1}{Q_{n}}\overset{n}{\underset{k=1}{\sum}}q_{n-k}S_{k}f,$

where

$Q_{n}:=\sum_{k=0}^{n-1}q_{k}.$

Here $\\{q_{k}:k\geq 0\\}$ is a sequence of nonnegative numbers, where
$q_{0}>0$ and

$\lim_{n\rightarrow\infty}Q_{n}=\infty.$

Then the summability method (9) generated by $\\{q_{k}:k\geq 0\\}$ is regular
if and only if (see [18])

$\underset{n\rightarrow\infty}{\lim}\frac{q_{n-1}}{Q_{n}}=0.$

In this paper we investigate regular Nörlund means only.

It is well-known (for details see e.g. [28]) that every Nörlund summability
method generated by non-increasing sequence $(q_{k},k\in\mathbb{N})$ is
regular, but Nörlund means generated by non-decreasing sequence
$(q_{k},k\in\mathbb{N})$ is not always regular.

The representation

$t_{n}f\left(x\right)=\underset{G_{m}}{\int}f\left(t\right)F_{n}\left(x-t\right)d\mu\left(t\right)$

play central roles in the sequel, where

(10) $F_{n}=:\frac{1}{Q_{n}}\overset{n}{\underset{k=1}{\sum}}q_{n-k}D_{k}$

is called the kernels of the Nörlund means.

If we invoke Abel transformation we get the following identities:

(11) $\displaystyle
Q_{n}:=\overset{n-1}{\underset{j=0}{\sum}}q_{j}=\overset{n}{\underset{j=1}{\sum}}q_{n-j}\cdot
1=\overset{n-1}{\underset{j=1}{\sum}}\left(q_{n-j}-q_{n-j-1}\right)j+q_{0}n$

and

(12)
$t_{n}f=\frac{1}{Q_{n}}\left(\overset{n-1}{\underset{j=1}{\sum}}\left(q_{n-j}-q_{n-j-1}\right)j\sigma_{j}f+q_{0}n\sigma_{n}f\right).$

##  3\. Formulation of Main Results

Based on estimate (1) we can prove our next main results:

######  Theorem 1.

Let $M_{N}\leq n<M_{N+1}$ and ${{t}_{n}}$ be a regular Nörlund mean generated
by non-decreasing sequence $\\{q_{k}:k\in\mathbb{N}\\},$ in sign
$q_{k}\uparrow.$ Then, for some $f\in L^{p}(G_{m}),$ where $1\leq p<\infty,$

$\|t_{n}f-f\|_{p}\leq\frac{3R^{3}}{Q_{n}}\sum_{i=0}^{N-1}M_{i}q_{n-M_{i}}\omega_{p}\left(\frac{1}{M_{i}},f\right)+2R^{3}\omega_{p}\left(\frac{1}{M_{N}},f\right).$

######  Theorem 2.

Let ${{t}_{n}}$ be Nörlund mean generated by non-increasing sequence
$\\{q_{k}:k\in\mathbb{N}\\}$, in sign $q_{k}\downarrow$. Then, for some $f\in
L^{p}(G_{m}),$ where $1\leq p<\infty,$

$\displaystyle\|t_{M_{n}}f-f\|_{p}\leq
3R^{2}\sum_{s=0}^{n}\frac{M_{s}}{M_{n}}\omega_{p}\left(1/M_{s},f\right)+C\overset{n-1}{\underset{s=0}{\sum}}\frac{(n-s)M_{s}}{M_{n}}\frac{q_{M_{s}}}{q_{M_{n}}}\omega_{p}\left(1/M_{s},f\right).$

######  Theorem 3.

Let $M_{N}\leq n<M_{N+1}$ and ${{t}_{n}}$ be Nörlund mean generated by non-
increasing sequence $\\{q_{k}:k\in\mathbb{N}\\}$, in sign $q_{k}\downarrow$,
satisfying the condition

(13) $\frac{1}{Q_{n}}=O\left(\frac{1}{n}\right),\text{\ \ as \ \
}n\rightarrow\infty.$

Then, for some $f\in L^{p}(G_{m}),$ where $1\leq p<\infty,$

$\displaystyle\|t_{n}f-f\|_{p}\leq
C\sum_{j=0}^{N}\frac{M_{j}}{M_{N}}\omega_{p}\left(1/M_{j},f\right).$

As a consequence we obtain the following similar result proved in Móricz and
Siddiqi [19]:

######  Corollary 1.

Let $\\{q_{k}:k\geq 0\\}$ be a sequence of non-negative numbers such that in
case $q_{k}\uparrow$ condition

(14) $\frac{q_{n-1}}{Q_{n}}=O\left(\frac{1}{n}\right),\text{\ \ as \ \
}n\rightarrow\infty.$

is satisfied, while in case $q_{k}\downarrow$ condition (13) is satisfied. If
$f\in Lip(\alpha,p)$ for some $\alpha>0$ and $1\leq p<\infty,$ then

(15) $\|t_{n}f-f\|_{p}=\left\\{\begin{array}[]{l}\text{
}O(n^{-\alpha}),\text{\qquad\qquad if\qquad}0<\alpha<1,\\\ \text{
}O(n^{-1}\log n),\text{\qquad if \quad\ }\alpha=1,\\\ \text{
}O(n^{-1}),\text{\qquad\qquad if \quad\ \ }\alpha>1,\end{array}\right.$

As a consequence we obtain the following similar result proved in Móricz and
Siddiqi [19]:

######  Corollary 2.

a) Let ${{t}_{n}}$ be Nörlund means generated by non-decreasing sequence
$\\{q_{k}:k\in\mathbb{N}\\}$ satisfying regularity condition (14). If $f\in
Lip(\alpha,p)$ for some $\alpha>0$ and $1\leq p<\infty,$ then ${t}_{n}f$
converge to $f$ in $L_{p}(G_{m})$ norm.

b) Let ${{t}_{n}}$ be Nörlund mean generated by non-increasing sequence
$\\{q_{k}:k\in\mathbb{N}\\}$ satisfying condition (13). If $f\in
Lip(\alpha,p)$ for some $\alpha>0$ and $1\leq p<\infty,$ then ${t}_{n}f$
converge to $f$ in $L_{p}(G_{m})$ norm.

##  4\. Proofs

###### Proof of Theorem 1..

Let $M_{N}\leq n<M_{N+1}.$ Since $t_{n}$ be regular Nörlund means generated by
the sequence of non-decreasing sequence $\\{q_{k}:k\in\mathbb{N}\\},$ by
combining (11) and (12) we can conclude that

$\displaystyle\|t_{n}f(x)-f(x)\|_{p}$ $\displaystyle\leq$
$\displaystyle\frac{1}{Q_{n}}\left(\overset{n-1}{\underset{j=1}{\sum}}\left(q_{n-j}-q_{n-j-1}\right)j\|\sigma_{j}f(x)-f(x)\|_{p}+q_{0}n\|\sigma_{n}f(x)-f(x)\|_{p}\right)$
$\displaystyle:=$ $\displaystyle I+II.$

Furthermore,

$\displaystyle I$ $\displaystyle=$
$\displaystyle\frac{1}{Q_{n}}\overset{M_{N}-1}{\underset{j=1}{\sum}}\left(q_{n-j}-q_{n-j-1}\right)j\|\sigma_{j}f(x)-f(x)\|_{p}$
$\displaystyle+$
$\displaystyle\frac{1}{Q_{n}}\overset{n-1}{\underset{j=M_{N}}{\sum}}\left(q_{n-j}-q_{n-j-1}\right)j\|\sigma_{j}f(x)-f(x)\|_{p}:=I_{1}+I_{2}.$

Now we estimate each terms separately. By applying (1) for $I_{1}$ we can
conclude that

$\displaystyle I_{1}$ $\displaystyle\leq$
$\displaystyle\frac{R^{2}}{Q_{n}}\overset{N-1}{\underset{k=0}{\sum}}\overset{M_{k+1}-1}{\underset{j=M_{k}}{\sum}}\left(q_{n-j}-q_{n-j-1}\right)j\sum_{s=0}^{k}\frac{M_{s}}{M_{k}}\omega_{p}\left(1/M_{s},f\right)$
$\displaystyle\leq$
$\displaystyle\frac{R^{2}}{Q_{n}}\overset{N-1}{\underset{k=0}{\sum}}M_{k+1}\overset{M_{k+1}-1}{\underset{j=M_{k}}{\sum}}\left(q_{n-j}-q_{n-j-1}\right)\sum_{s=0}^{k}\frac{M_{s}}{M_{k}}\omega_{p}\left(1/M_{s},f\right)$
$\displaystyle\leq$
$\displaystyle\frac{R^{3}}{Q_{n}}\overset{N-1}{\underset{k=0}{\sum}}\left(q_{n-M_{k}}-q_{n-M_{k+1}}\right)\sum_{s=0}^{k}M_{s}\omega_{p}\left(1/M_{s},f\right)$
$\displaystyle\leq$
$\displaystyle\frac{R^{3}}{Q_{n}}\overset{N-1}{\underset{s=0}{\sum}}M_{s}\omega_{p}\left(1/M_{s},f\right)\sum_{k=s}^{N-1}\left(q_{n-M_{k}}-q_{n-M_{k+1}}\right)$
$\displaystyle\leq$
$\displaystyle\frac{R^{3}}{Q_{n}}\overset{N-1}{\underset{s=0}{\sum}}M_{s}q_{n-M_{s}}\omega_{p}\left(1/M_{s},f\right).$

It is evident that

$\displaystyle\quad\quad I_{2}$ $\displaystyle\leq$
$\displaystyle\frac{R^{2}}{Q_{n}}\overset{n-1}{\underset{j=M_{N}}{\sum}}\left(q_{n-j}-q_{n-j-1}\right)j\sum_{s=0}^{N}\frac{M_{s}}{M_{N}}\omega_{p}\left(1/M_{s},f\right)$
$\displaystyle\leq$
$\displaystyle\frac{R^{2}M_{N+1}}{Q_{n}}\overset{n-1}{\underset{j=M_{N}}{\sum}}\left(q_{n-j}-q_{n-j-1}\right)\sum_{s=0}^{N}\frac{M_{s}}{M_{N}}\omega_{p}\left(1/M_{s},f\right)$
$\displaystyle\leq$
$\displaystyle\frac{R^{3}q_{n-M_{N}}}{Q_{n}}\sum_{s=0}^{N}M_{s}\omega_{p}\left(1/M_{s},f\right)$
$\displaystyle\leq$
$\displaystyle\frac{R^{3}}{Q_{n}}\overset{N}{\underset{s=0}{\sum}}M_{s}q_{n-M_{s}}\omega_{p}\left(1/M_{s},f\right).$

For $II$ we have that

$\displaystyle II$ $\displaystyle\leq$
$\displaystyle\frac{q_{0}R^{2}M_{N+1}}{Q_{n}}\sum_{s=0}^{N}\frac{M_{s}}{M_{N}}\omega_{p}\left(1/M_{s},f\right)$
$\displaystyle\leq$
$\displaystyle\frac{R^{3}}{Q_{n}}\sum_{s=0}^{N-1}M_{s}q_{n-M_{s}}\omega_{p}\left(1/M_{s},f\right)+R^{3}\omega_{p}\left(1/M_{N},f\right).$

The proof is complete. ∎

###### Proof of Theorem 2..

By using (2) we find that

(18) $\displaystyle t_{M_{n}}f=D_{M_{n}}\ast
f-\frac{1}{Q_{M_{n}}}\overset{M_{n}-1}{\underset{k=0}{\sum}}q_{k}\left(\left(\psi_{M_{n}-1}\overline{D_{k}}\right)\ast
f\right).$

By using Abel transformation we get

$\displaystyle t_{M_{n}}f$ $\displaystyle=$ $\displaystyle D_{M_{n}}\ast
f-\frac{1}{Q_{M_{n}}}\overset{M_{n}-2}{\underset{j=0}{\sum}}\left(q_{j}-q_{j+1}\right)j((\psi_{M_{n}-1}\overline{K_{j}})\ast
f)$ $\displaystyle-$
$\displaystyle\frac{1}{Q_{M_{n}}}q_{M_{n}-1}(M_{n}-1)(\psi_{M_{n}-1}\overline{K}_{M_{n}-1}\ast
f)$ $\displaystyle=$ $\displaystyle D_{M_{n}}\ast
f-\frac{1}{Q_{M_{n}}}\overset{M_{n}-2}{\underset{j=0}{\sum}}\left(q_{j}-q_{j+1}\right)j((\psi_{M_{n}-1}\overline{K_{j}})\ast
f)$ $\displaystyle-$
$\displaystyle\frac{1}{Q_{M_{n}}}q_{M_{n}-1}M_{n}(\psi_{M_{n}-1}\overline{K}_{M_{n}}\ast
f)$ $\displaystyle+$
$\displaystyle\frac{q_{M_{n}-1}}{Q_{M_{n}}}(\psi_{M_{n}-1}\overline{D}_{M_{n}}\ast
f)$

and

$\displaystyle t_{M_{n}}f(x)-f(x)$ $\displaystyle=$
$\displaystyle\int_{G_{m}}(f(x+t)-f(x))D_{M_{n}}(t)dt$ $\displaystyle-$
$\displaystyle\frac{1}{Q_{M_{n}}}\overset{M_{n}-2}{\underset{j=0}{\sum}}\left(q_{j}-q_{j+1}\right)j\int_{G_{m}}\left(f(x+t)-f(x)\right)\psi_{M_{n}-1}(t)\overline{K}_{j}(t)dt$
$\displaystyle-$
$\displaystyle\frac{1}{Q_{M_{n}}}q_{M_{n}-1}M_{n}\int_{G_{m}}(f(x+t)-f(x))\psi_{M_{n}-1}(t)\overline{K}_{M_{n}}(t)dt$
$\displaystyle+$
$\displaystyle\frac{q_{M_{n}-1}}{Q_{M_{n}}}\int_{G_{m}}(f(x+t)-f(x))\psi_{M_{n}-1}(t)\overline{D}_{M_{n}}(t)dt$
$\displaystyle=$ $\displaystyle I+II+III+IV.$

By combining generalized Minkowski’s inequality and (4) we find that

$\|I\|_{p}\leq\int_{I_{n}}\|f(x+t)-f(x))\|_{p}D_{M_{n}}(t)dt\leq\omega_{p}\left(1/M_{n},f\right).$

and

$\|IV\|_{p}\leq\int_{I_{n}}\|f(x+t)-f(x))\|_{p}{D}_{M_{n}}(t)dt\leq\omega_{p}\left(1/M_{n},f\right).$

Since $M_{n}q_{M_{n}-1}\leq Q_{M_{n}},$ for any $n\in\mathbb{N},$ if combine
(8) and generalized Minkowski’s inequality we get

$\displaystyle\|III\|_{p}$ $\displaystyle\leq$
$\displaystyle\int_{G_{m}}\left\|f\left(x+t\right)-f\left(x\right)\right\|_{p}\left|\overline{K}_{M_{n}}\left(t\right)\right|d\mu(t)$
$\displaystyle=$
$\displaystyle\int_{I_{n}}\left\|f\left(x+t\right)-f\left(x\right)\right\|_{p}\left|\overline{K}_{M_{n}}\left(t\right)\right|d\mu(t)$
$\displaystyle+$
$\displaystyle\sum_{s=0}^{n-1}\sum_{n_{s}=1}^{m_{s}-1}\int_{I_{n}\left(n_{s}e_{s}\right)}\left\|f\left(x+t\right)-f\left(x\right)\right\|_{p}\left|\overline{K}_{M_{n}}\left(t\right)\right|d\mu(t)$
$\displaystyle\leq$
$\displaystyle\int_{I_{n}}\left\|f\left(x+t\right)-f\left(x\right)\right\|_{p}\frac{M_{n}+1}{2}d\mu(t)$
$\displaystyle+$
$\displaystyle\sum_{s=0}^{n-1}M_{s+1}\sum_{n_{s}=1}^{m_{s}-1}\int_{I_{n}\left(n_{s}e_{s}\right)}\left\|f\left(x+t\right)-f\left(x\right)\right\|_{p}d\mu(t)$
$\displaystyle\leq$
$\displaystyle\omega_{p}\left(1/M_{n},f\right)\int_{I_{n}}\frac{M_{n}+1}{2}d\mu(t)$
$\displaystyle+$
$\displaystyle\sum_{s=0}^{n-1}M_{s+1}\sum_{n_{s}=1}^{m_{s}-1}\int_{I_{n}\left(n_{s}e_{s}\right)}\omega_{p}\left(1/M_{s},f\right)d\mu(t)$
$\displaystyle\leq$
$\displaystyle\omega_{p}\left(1/M_{n},f\right)+R^{2}\sum_{s=0}^{n-1}\frac{M_{s}}{M_{n}}\omega_{p}\left(1/M_{s},f\right)\leq
R^{2}\sum_{s=0}^{n}\frac{M_{s}}{M_{n}}\omega_{p}\left(1/M_{s},f\right).$

This estimate also follows that

(21)
$M_{n}\int_{G_{m}}\left\|f\left(x+t\right)-f\left(x\right)\right\|_{p}|\overline{K}_{M_{n}}(t)|\,d\mu(t)\leq
R^{2}\sum_{s=0}^{n}M_{s}\omega_{p}\left(1/M_{s},f\right).$

Let $M_{k}\leq j<M_{k+1}$ By applying (6) and (21) we find that

$\displaystyle
j\int_{G_{m}}\left\|f\left(x+t\right)-f\left(x\right)\right\|_{p}|\overline{K}_{j}(t)|d\mu(t)$
$\displaystyle\leq$ $\displaystyle
C\sum_{l=0}^{k}\sum_{s=0}^{l}M_{s}\omega_{p}\left(1/M_{s},f\right).$

Hence, by combining (6) and (4) we find that

$\displaystyle\|II\|_{p}$ $\displaystyle\leq$
$\displaystyle\frac{1}{Q_{M_{n}}}\overset{M_{n}-1}{\underset{j=0}{\sum}}\left(q_{j}-q_{j+1}\right)j\int_{G_{m}}\|f(x+t)-f(x)\|_{p}|\overline{K}_{j}(t)|d\mu(t)$
$\displaystyle\leq$
$\displaystyle\frac{1}{Q_{M_{n}}}\overset{n-1}{\underset{k=0}{\sum}}\overset{M_{k+1}-1}{\underset{j=M_{k}}{\sum}}\left(q_{j}-q_{j+1}\right)j\int_{G_{m}}\|f(x+t)-f(x)\|_{p}|\overline{K}_{j}(t)|d\mu(t)$
$\displaystyle\leq$
$\displaystyle\frac{C}{Q_{M_{n}}}\overset{n-1}{\underset{k=0}{\sum}}\overset{M_{k+1}-1}{\underset{j=M_{k}}{\sum}}\left(q_{j}-q_{j+1}\right)\sum_{l=0}^{k}\sum_{s=0}^{l}M_{s}\omega_{p}\left(1/M_{s},f\right)$
$\displaystyle\leq$
$\displaystyle\frac{C}{Q_{M_{n}}}\overset{n-1}{\underset{k=0}{\sum}}\left(q_{M_{k}}-q_{M_{k+1}}\right)\sum_{l=0}^{k}\sum_{s=0}^{l}M_{s}\omega_{p}\left(1/M_{s},f\right)$
$\displaystyle\leq$
$\displaystyle\frac{C}{Q_{M_{n}}}\overset{n-1}{\underset{l=0}{\sum}}\sum_{k=l}^{n-1}\left(q_{M_{k}}-q_{M_{k+1}}\right)\sum_{s=0}^{l}M_{s}\omega_{p}\left(1/M_{s},f\right)$
$\displaystyle\leq$
$\displaystyle\frac{C}{Q_{M_{n}}}\overset{n-1}{\underset{l=0}{\sum}}q_{M_{l}}\sum_{s=0}^{l}M_{s}\omega_{p}\left(1/M_{s},f\right)$
$\displaystyle\leq$
$\displaystyle\frac{C}{Q_{M_{n}}}\overset{n-1}{\underset{s=0}{\sum}}M_{s}\omega_{p}\left(1/M_{s},f\right)\sum_{l=s}^{n-1}q_{M_{l}}$
$\displaystyle\leq$
$\displaystyle\frac{C}{Q_{M_{n}}}\overset{n-1}{\underset{s=0}{\sum}}M_{s}\omega_{p}\left(1/M_{s},f\right)q_{M_{s}}(n-s)$
$\displaystyle\leq$ $\displaystyle
C\overset{n-1}{\underset{s=0}{\sum}}\frac{(n-s)M_{s}}{M_{n}}\frac{q_{M_{s}}}{q_{M_{n}}}\omega_{p}\left(1/M_{s},f\right).$

The proof is complete. ∎

###### Proof of theorem 3..

Let $M_{N}\leq n<M_{N+1}.$ Since $t_{n}$ be regular Nörlund means, generated
by sequence of non-increasing numbers $\\{q_{k}:k\in\mathbb{N}\\}$ by
combining (11) and (12), we can conclude that

$\displaystyle\|t_{n}f(x)-f(x)\|_{p}$ $\displaystyle\leq$
$\displaystyle\frac{1}{Q_{n}}\left(\overset{n-1}{\underset{j=1}{\sum}}\left(q_{n-j-1}-q_{n-j}\right)j\|\sigma_{j}f(x)-f(x)\|_{p}+q_{0}n\|\sigma_{n}f(x)-f(x)\|_{p}\right)$
$\displaystyle:=$ $\displaystyle I+II.$

Furthermore,

$\displaystyle I$ $\displaystyle=$
$\displaystyle\frac{1}{Q_{n}}\overset{M_{N}-1}{\underset{j=1}{\sum}}\left(q_{n-j-1}-q_{n-j}\right)j\|\sigma_{j}f(x)-f(x)\|_{p}$
$\displaystyle+$
$\displaystyle\frac{1}{Q_{n}}\overset{n-1}{\underset{j=M_{N}}{\sum}}\left(q_{n-j-1}-q_{n-j}\right)j\|\sigma_{j}f(x)-f(x)\|_{p}$
$\displaystyle=$ $\displaystyle I_{1}+I_{2}.$

Analogously to (4) we get that

$\displaystyle I_{1}$ $\displaystyle\leq$
$\displaystyle\frac{R^{3}}{Q_{n}}\overset{N-1}{\underset{k=0}{\sum}}\left(q_{n-M_{k+1}}-q_{n-M_{k}}\right)\sum_{s=0}^{k}M_{s}\omega_{p}\left(1/M_{s},f\right)$
$\displaystyle\leq$
$\displaystyle\frac{R^{3}}{Q_{n}}\sum_{s=0}^{N-1}M_{s}\omega_{p}\left(1/M_{s},f\right)\overset{N-1}{\underset{k=s}{\sum}}\left(q_{n-M_{k+1}}-q_{n-M_{k}}\right)$
$\displaystyle=$
$\displaystyle\frac{R^{3}}{Q_{n}}\sum_{s=0}^{N-1}M_{s}\omega_{p}\left(1/M_{s},f\right)(q_{n-M_{N}}-q_{n-M_{s}})$
$\displaystyle\leq$
$\displaystyle\frac{R^{3}q_{n-M_{N}}}{Q_{n}}\sum_{s=0}^{N-1}M_{s}\omega_{p}\left(1/M_{s},f\right)\leq\frac{R^{3}q_{0}}{Q_{n}}\sum_{s=0}^{N-1}M_{s}\omega_{p}\left(1/M_{s},f\right).$

Analogously to (4) we find that

$\displaystyle I_{2}$ $\displaystyle\leq$
$\displaystyle\frac{R^{2}}{Q_{n}}\overset{n-1}{\underset{j=1}{\sum}}\left(q_{n-j-1}-q_{n-j}\right)j\sum_{s=0}^{N}\frac{M_{s}}{M_{N}}\omega_{p}\left(1/M_{s},f\right)$
$\displaystyle=$
$\displaystyle\frac{R^{2}}{Q_{n}}\left(nq_{0}-Q_{n}\right)\sum_{s=0}^{N}\frac{M_{s}}{M_{N}}\omega_{p}\left(1/M_{s},f\right)$
$\displaystyle\leq$
$\displaystyle\frac{M_{N+1}R^{2}q_{0}}{Q_{n}M_{N}}\sum_{s=0}^{N}M_{s}\omega_{p}\left(1/M_{s},f\right)$
$\displaystyle\leq$
$\displaystyle\frac{R^{3}q_{0}}{Q_{n}}\sum_{s=0}^{N}M_{s}\omega_{p}\left(1/M_{s},f\right).$

For $II$ we have that

$\displaystyle II$ $\displaystyle\leq$
$\displaystyle\frac{q_{0}R^{2}M_{N+1}}{Q_{n}}\sum_{s=0}^{N}\frac{M_{s}}{M_{N}}\omega_{p}\left(1/M_{s},f\right)$
$\displaystyle\leq$
$\displaystyle\frac{R^{3}q_{0}}{Q_{n}}\sum_{s=0}^{N}M_{s}\omega_{p}\left(1/M_{s},f\right).$

Using (13) we obtain estimate above so the proof is complete. ∎

## References

  * [1] G. N. Agaev, N. Ya. Vilenkin, G. M. Dzhafarly and A. I. Rubinshtein, Multiplicative systems of functions and harmonic analysis on zero-dimensional groups, Baku, Ehim, 1981 (in Russian).
  * [2] N. Areshidze and G. Tephnadze, Approximation by Nörlund means with respect to Walsh system in Lebesgue spaces, Math. Inequal. Appl., (to appear).
  * [3] D. Baramidze, N. Nadirashvili, L.-E. Persson and G. Tephnadze, Some weak-type inequalities and almost everywhere convergence of Vilenkin-Nörlund means, J. Inequal. Appl., 2023, 2023 (1), 66.
  * [4] L. Baramidze, L. E. Persson, G. Tephnadze and P. Wall, Sharp $H_{p}-L_{p}$ type inequalities of weighted maximal operators of Vilenkin-Nörlund means and its applications, J. Inequal. Appl., 2016, DOI: 10.1186/s13660-016-1182-1, 2016, 2016(1), 242.
  * [5] I. Blahota, K. Nagy, G. Tephnadze, Approximation by Marcinkiewicz $\Theta$-means of double Walsh-Fourier series, Math. Inequal. Appl., 22, 3 (2019) 837-853.
  * [6] I. Blahota, G. Tephnadze, On the $(C,\alpha)$-means with respect to the Walsh system, Anal. Math., 40 (2014), 161-174.
  * [7] I. Blahota, G. Tephnadze, R. Toledo, Strong convergence theorem of $(C,\alpha)$-means with respect to the Walsh system, Tohoku Math. J., 67, 4 (2015), 573-584.
  * [8] I. Blahota, L.E. Persson, G. Tephnadze, On the Nörlund means of Vilenkin-Fourier series, Czech. Math J., 65, 4 (2015), 983-1002.
  * [9] I. Blahota and K. Nagy, Approximation by $\Theta$-means of Walsh-Fourier series, Anal. Math. 44, 1 (2018), 57-71.
  * [10] I. Blahota, K. Nagy and G. Tephnadze, Approximation by Marcinkiewicz $\Theta$-means of double Walsh-Fourier series, Math. Inequal. Appl., 22, 3 (2019), 837-853.
  * [11] S. Fridli, On the rate of convergence of Cesaro means of Walsh-Fourier series, J. Approx. Theory, 76, 1 (1994), 31-53.
  * [12] S. Fridli, P, Manchanda and A. Siddiqi, Approximation by Walsh-Nörlund means, Acta Sci. Math., 74 (3-4), 593-608.
  * [13] G. Gát, Cesàro means of integrable functions with respect to unbounded Vilenkin systems. J. Approx. Theory, 124, 1 (2003), 25-43.
  * [14] U. Goginava, Maximal operators of Walsh-Nörlund means on the dyadic Hardy spaces, Acta Math, Hungar., 169, 1 (2023), 171-190.
  * [15] U. Goginava, Marcinkiewicz-Fejer means of d-dimensional Walsh-Fourier series, 307, 1 (2015), 206 - 218.
  * [16] B. I. Golubov, A. V. Efimov and V. A. Skvortsov, Walsh series and transforms, Kluwer Academic Publishers Group, Dordrecht, 1991.
  * [17] N. Memić, An Estimate of the Maximal Operator of the Nörlund Logarithmic Means with Respect to the Character System of the Group of 2-adic Integers on the Hardy Space $H_{1},$ 48, 6 (2022), 3381-3391.
  * [18] C. N. Moore, Summable series and convergence factors, Summable series and convergence factors. Dover Publications, Inc., New York, 1966\.
  * [19] F. Móricz and A. Siddiqi, Approximation by Nörlund means of Walsh-Fourier series, J. Approx. Theory, 70, 3 (1992), 375-389.
  * [20] K. Nagy, Approximation by Nörlund means of double Walsh-Fourier series for Lipschitz functions, Math. Inequal. Appl., 15, 2 (2012), 301-322.
  * [21] K. Nagy, Approximation by Nörlund means of Walsh-Kaczmarz-Fourier series. Georgian Math. J. 18, 1 (2011), 147-162.
  * [22] K. Nagy, Approximation by Nörlund means of quadratical partial sums of double Walsh-Fourier series, Anal. Math., 36, 4 (2010), 299-319.
  * [23] K. Nagy and G. Tephnadze, Walsh-Marcinkiewicz means and Hardy spaces, Cent. Eur. J. Math., 12, 8 (2014), 1214-1228.
  * [24] K. Nagy and G. Tephnadze, Approximation by Walsh-Marcinkiewicz means on the Hardy space , Kyoto J. Math., 54, 3 (2014), 641-652.
  * [25] K. Nagy and G. Tephnadze, Strong convergence theorem for Walsh-Marcinkiewicz means, Math. Inequal. Appl., 19, 1 (2016), 185-195.
  * [26] K. Nagy and G. Tephnadze, The Walsh-Kaczmarz-Marcinkiewicz means and Hardy spaces, Acta math. Hung., 149, 2 (2016), 346-374.
  * [27] J. Pál and P. Simon, On a generalization of the concept of derivative, Acta Math. Hung., 29 (1977), 155-164.
  * [28] L. E. Persson, G. Tephnadze and F. Weisz, Martingale Hardy Spaces and Summability of Vilenkin-Fourier Series, book manuscript, Springer, 2022.
  * [29] F. Schipp, Certain rearranngements of series in the Walsh series, Mat. Zametki, 18 (1975), 193-201.
  * [30] F. Schipp, W. R. Wade, P. Simon and J. Pál, Walsh series. An introduction to dyadic harmonic analysis, Adam Hilger, Ltd., Bristol, 1990.
  * [31] P Simon, Strong convergence of certain means with respect to the Walsh-Fourier series, Acta Math. Hungar., 49, 3-4 (1987), 425-431.
  * [32] P Simon, On the Cesáro summability with respect to the Walsh?Kaczmarz system, J. Approx. Theory, 106, 2 (2000), 249-261.
  * [33] G. Tutberidze, Modulus of continuity and boundedness of subsequences of Vilenkin-Fejer means in the martingale Hardy spaces, Geo. Math. J., 29, 1 (2022), 153-162.
  * [34] N. Y. Vilenkin, On a class of complete orthonormal systems, Amer. Math. Soc. Transi., 28, (2) (1963), 1-35.
  * [35] N. Y. Vilenkin, On the theory of lacunary orthogonal systems, (Russian) Izv. Akad. Nauk SSSR, Ser. Mat., 13 (1949), 245-252.
  * [36] N. Y. Vilenkin, On the theory of Fourier integrals on topological groups, Mat. Sbornik N. S., 30, (72) (1952), 233-244.
  * [37] F. Weisz, Martingale Hardy Spaces and their Applications in Fourier Analysis, Springer, Berlin-Heideiberg-New York, 1994.
  * [38] F. Weisz, Hardy spaces and Cesàro means of two-dimensional Fourier series, Bolyai Soc. Math. Studies, (1996), 353-367.
  * [39] A. Zygmund, Trigonometric Series, Vol. 1, Cambridge Univ. Press, 1959.