# Encrypted Simultaneous Control of Joint Angle and Stiffness of Antagonistic
Pneumatic Artificial Muscle Actuator by Polynomial Approximation

Yuta Takeda, Takaya Shin, Kaoru Teranishi, and Kiminao Kogiso This work was
supported by JSPS KAKENHI Grant Numbers JP22H01509 and JP21K19762.Y. Takeda,
K. Teranishi, and K. Kogiso are with the Department of Mechanical and
Intelligent Systems Engineering, The University of Electro-Communications,
1-5-1 Chofugaoka, Chofu, Tokyo 1828585, Japan. e-mail<EMAIL_ADDRESS>Shin is with DAIHEN, Inc.

###### Abstract

This study proposes an encrypted simultaneous control system for an
antagonistic pneumatic artificial muscle (PAM) actuator toward developing a
cybersecure and flexible actuator. First, a novel simultaneous control system
design is considered for the joint angle and stiffness of a PAM actuator in a
model-based design approach, facilitating the use of an encrypted control
method. The designed controller includes a contraction force model expressed
as rational polynomial functions, which makes it difficult to encrypt the
controller. To overcome this difficulty, a least absolute shrinkage and
selection operator (LASSO)-based polynomial approximation is employed for a
rational controller. The resulting polynomial controller is then transformed
into a matrix-vector product form, which enables the use of a specific
homomorphic encryption scheme to develop an encrypted simultaneous control
system for the PAM actuator. Finally, this study quantitatively evaluates the
tracking control performance of the original, approximated, and encrypted
controllers. The experimental results show that the proposed encrypted
controller achieves simultaneous tracking of the joint angle and stiffness
with a tracking error of less than 2.7 %.

###### Index Terms:

PAM, encrypted control, polynomial-type controller, LASSO, simultaneous
control, experimental validation

##  I Introduction

The McKibben pneumatic artificial muscle (PAM) was first developed in the
1950s [1]. A PAM is inflated by injecting compressed air into a rubber tube
wrapped in a nonstretchable mesh to generate a contraction force. This
structure offers several advantages over traditional motors and cylinders,
including compact size, lightweight, and flexibility [2, 3]. Because PAMs
produce contraction force in a single direction, an antagonistic configuration
is typically adopted for practical applications, where two PAMs are arranged
in parallel through joints [4, 5, 6, 7, 8, 9, 10, 11, 12]. This structure
facilitates rotational motion and closely mimics human movements [3]. Owing to
these benefits, a review paper [13] stated that PAMs have found considerable
use in the development of rehabilitation and assistive devices that can be
comfortably worn and operated by users. Furthermore, numerous studies have
explored stiffness or compliance control [14, 15, 16, 17, 18] and developed
methods of simultaneous control of the joint angle and stiffness to leverage
the potential of PAMs fully [19, 20, 21, 22, 23, 24]. These active controls
enhance safety during contact and collision between robots and humans or their
environments while also augmenting the comfort of wearable devices.

Applications that integrate wireless network technology into PAM actuator
systems include remote surgical robots [25], remotely operated cranes [26],
teleoperation of pneumatic robots [27], and walking-assistive devices [28,
29]. The use of wireless communications enhances the portability of devices
and promotes the self-rehabilitation of patients at home. In self-
rehabilitation, the monitoring of wearable devices through a network by
physical therapists helps prevent injury due to falls. Meanwhile, attention
must be paid to the cybersecurity of networked actuator systems. Cyberattack
incidents and risk analyses have been reported, such as the falsification of
control parameters by Stuxnet to destroy process plants [30], an unmanned
aerial vehicle compromised by hijacking video streaming [31], and the
compromising of a robot controller [32].

To develop secure PAM actuator systems, networked actuator systems must be
equipped with cybersecurity countermeasures. Encrypted control was proposed as
a cyberserucrity countermeasure for networked control systems in [33] by
integrating the homomorphism of a specific public-key encryption scheme into a
linear or polynomial-type control system. Encrypted control conceals the
control parameters inside the control device and the signals over the
communication links, resulting in protection against eavesdropping [33, 34,
35, 36] and real-time detection of cyberattacks [37, 33, 38]. Encrypted
implementation is expected to enhance the cybersecurity of control systems.

However, designing a nonlinear controller that is not a polynomial, such as a
rational polynomial or switching controllers, is difficult. A simultaneous
control system for the joint angle and stiffness of an antagonistic PAM
actuator uses a rational polynomial-type controller [19, 22, 23]. In
additional, encrypted proportional-integral (PI) control of the PAM joint
angle or torque was considered in [36], and in [39], the implementation of
encrypted polynomial controllers was demonstrated, verifying its effectiveness
through a numerical toy example. Unfortunately, to the best of our knowledge,
no studies have been conducted on encrypting such nonlinear controllers.

The objective of this study is to propose an encrypted simultaneous control
system for an antagonistic PAM actuator to develop a cybersecure and flexible
actuator. In this study, the simultaneous tracking control of the joint angle
and stiffness of the antagonistic PAM actuator is considered for a step-like
reference. We present a novel model-based nonlinear control system comprising
three components: a reference generator, a contraction force estimator, and PI
controllers. However, a model-based controller includes specific rational
functions, and applying the encrypted control method is challenging. To
overcome this difficulty, this study approximates the designed controller as a
polynomial-type controller that is friendly to controller encryption. We then
quantitatively investigate the effects of polynomial approximation and
controller encryption on the control performance to evaluate the developed
encrypted control system. Finally, an experimental investigation confirms that
the proposed encrypted control system enables simultaneous tracking to the
step-like reference within an acceptable tracking error.

This study makes the following main contributions: it is the first to present
an encrypted nonlinear control system for simultaneously tracking the joint
angle and stiffness of an antagonistic PAM actuator, and it provides a novel
form of the contraction force model, which helps reduce the number of
approximations that cause control performance degradation throughout the
controller encryption procedure.

The remainder of this paper is organized as follows. Section II introduces the
antagonistic PAM actuator with a control objective and presents a model-based
simultaneous control method for tracking the joint angle and stiffness of the
PAM actuator to the reference. Section III discusses the polynomial
approximation to acquire a friendly form of the controller to realize the
encrypted control. Section IV proposes an encrypted-controlled PAM actuator to
maintain the controller and communication secrets and investigates the impact
of the secure implementation on the control performance to highlight the
effectiveness of the proposed encrypted controller. Finally, Section V
concludes the study.

(a)

(b)

Figure 1: Antagonistic PAM actuator system considered [40]. (a) Schematic of the entire system and (b) geometric structure of the PAM joint [24]. TABLE I: Experimental equipments Name | Specifications
---|---
PAM | AirMuscle, Kanda Tsushin Kogyo,
| Length: 170 mm (be given a 2 kgf mass),
| Diameter: 0.5 inches
Valve (PDCV) | 5/3-way valve, FESTO,
| a critical frequency of 125 Hz.
Torque sensor | UTM II-10 Nm(R), UNIPULSE,
| range: $\pm 10$ Nm.
Encoder | UTM II-10 Nm(R), UNIPULSE,
| Incremental, resolution: 2000 P/R
Pressure sensor | E8F2-B10C, OMRON,
| Range: $0$–$1$ MPa.
Air compressor | 6-25, JUN-AIR,
| Tank: 25 L; displacement: 60 L/min.
Pressure tank | AST-25G, EARTH MAN, 25 L.
PC | Ubuntu12.04, Xenomai2.6.2.1 Patch,
| CPU: 3.2 GHz, memory: 8 GB.

##  II Nonlinear Control System Design for Antagonistic PAM Actuator

This section introduces the antagonistic PAM actuator, its simultaneous
control system desingn, and the numerical evaluation of the control
performance of the designed control system.

###  II-A Antagonistic PAM Actuator System

Fig. 1LABEL:sub@sfig:eqatg shows an antagonistic PAM system. An antagonistic
PAM actuator is a joint actuator driven by two PAMs. One side of each PAM is
connected by a joint. The system consists of two PAMs, two proportional
directional control valves (PDCVs), an air tank, pressure sensors, a torque
meter, a rotary encoder, and a control PC. The tank stores compressed air
connected to the PDCVs and PAMs using air tubes. The airflow regulated by the
PDCVs drives the PAMs, and the joint rotates. The system inputs are the
voltage commands to the two PDCVs ($u_{1}$ and $u_{2}$), and the measured
values are the joint angle $\theta$, torque $\tau$, and inner pressure of the
PAMs ($P_{1}$ and $P_{2}$). The range of the rotation angle is $\pm$ 25 deg,
and the range of the output torque is $\pm$ 3.0 Nm. The sampling period
$T_{\rm s}$ was set to 20 ms. Table I lists the details of the experimental
equipment used.

The control objective considered in this study is to track the joint angle and
stiffness of the PAM actuator simultaneously with respect to a given
reference. The goal of this study is to develop an encrypted simultaneous
tracking control system for the joint angle and stiffness of a PAM actuator to
mitigate the control performance degradation caused by secure implementation.

###  II-B Joint Stiffness

Joint stiffness is an index of the difficulty of a joint rotating against
external forces [23, 24]. This study employs the model of joint stiffness
presented in [24], following the introduction of the notations of the physical
variables regarding the PAM. The geometric relationships of the antagonistic
PAM actuator are illustrated in Fig. 1LABEL:sub@fig:apam, the lengths of the
two PAMs, $l_{1}$ and $l_{2}$, are respectively given by
$l_{1}(k)=L_{0}-\Delta L(k)$ and $l_{2}(k)=L_{0}+\Delta L(k)$, where
$k\in\mathbb{Z}^{+}:=\\{1,2,\cdots\\}$ is a step, $\Delta L(k)\approx
r\sin\theta(k)$ is the horizontal displacement of the two PAMs, $r$ is the
radius of the joint, and $L_{0}$ is the PAM length when the joint is at the
horizontal position. The contraction force of the PAM $F_{i}$, $\forall
i\in{\mathcal{I}}:=\\{1,2\\}$ can be expressed as a function of the inner
pressure and length as follows:

$\displaystyle
F_{i}(l_{i}(k),P_{i}(k))=a_{i}(l_{i}(k))P_{i}(k)+b_{i}(l_{i}(k)),\ \,\forall
i\in\mathcal{I},$ (1)

with $a_{i}(l)=p^{a_{1}}_{i}l+p^{a_{2}}_{i}$ and
$b_{i}(l)=p^{b_{1}}_{i}l+p^{b_{2}}_{i},\,\forall i\in\mathcal{I}$.

When the joint torque $\tau$ generated by the two PAMs is given by
$\tau(k)=r\cos\theta(k)\left(F_{1}(k)-F_{2}(k)\right)$, the joint stiffness
$K_{P}$ is defined as the partial differentiation of joint torque $\tau$ by
joint angle $\theta$ as follows [24]:

$\displaystyle K_{P}(k)$ $\displaystyle=r\sin\theta(k)(F_{1}(k)-F_{2}(k))$ (2)
$\displaystyle\quad+r^{2}\cos^{2}\theta(k)\left(\frac{F_{1}(k)-\alpha_{1}(k)}{l_{1}(k)}+\frac{F_{2}(k)-\alpha_{2}(k)}{l_{2}(k)}\right),$

where $\alpha_{i}(k):=p^{a_{2}}_{i}P_{i}(k)+p^{b_{2}}_{i},\,\forall
i\in\mathcal{I}$.

Figure 2: Block diagram of angle-stiffness control system.

###  II-C Angle-Stiffness Controller

This study proposes a novel nonlinear controller for the simultaneous control
of joint angle and stiffness, which is highlighted in red in Fig. 2, as
follows:

$\displaystyle x(k+1)$ $\displaystyle=A(k;\theta)x(k)+g(k;\zeta),$ (3a)
$\displaystyle\begin{bmatrix}u_{1}(k)\\\ u_{2}(k)\end{bmatrix}$
$\displaystyle=C(k;\theta)x(k)+h(k;\zeta)+\begin{bmatrix}\beta_{1}\\\
\beta_{2}\end{bmatrix},$ (3b)

where $x\in\mathbb{R}^{3}$ is the state, $u_{1}$ and $u_{2}$ are the outputs,
and $\zeta\in\mathbb{R}^{5}$ is the input, denoted by $\zeta:=[\,P_{1}\
\,P_{2}\ \,\theta\ \,\bar{\theta}\ \,\bar{K}_{P}\,]^{\mathrm{T}}$, where
$P_{1}$, $P_{2}$, and $\theta$ can be measured by the sensor, and
$\bar{\theta}$ and $\bar{K}_{P}$ are the references. $\beta_{1}$ and
$\beta_{2}$ are the bias voltages of each PDCV.
$A(\theta)\in\mathbb{R}^{3\times 3}$ and $C(\theta)\in\mathbb{R}^{2\times 3}$
are time-varying, $\theta$-dependent coefficients, and
$g:\mathbb{R}^{5}\rightarrow\mathbb{R}^{3}$ and
$f:\mathbb{R}^{5}\rightarrow\mathbb{R}^{3}$ are nonlinear functions of
$\zeta$. The coefficients and nonlinear functions are clarified after
introducing a reference generator, pressure into the contraction force, and PI
controllers, which are components of the proposed controller and are explained
as follows.

####  II-C1 Reference Generator

The reference generator outputs the reference signals of the contraction force
for two PAMs, denoted as $\bar{F}_{1}$ and $\bar{F}_{2}$, taking the measured
joint angle $\theta$, computed torque command $\tau_{c}$, and the given
reference of joint stiffness $\bar{K}_{P}$ as inputs. Using [24],
$\bar{F}_{1}$ and $\bar{F}_{2}$ are expressed as rational polynomial functions
of the inputs $\theta$, $\tau_{c}$, and $\bar{K}_{P}$:

$\displaystyle\bar{F}_{1}(k)=$ $\displaystyle\
\frac{1}{r^{2}\cos^{2}\theta(k)}\frac{l_{1}(k)l_{2}(k)}{l_{1}(k)+l_{2}(k)}\biggl{[}\bar{K}_{P}(k)\biggr{.}+\Bigl{(}\frac{r\cos\theta(k)}{l_{2}(k)}-$
$\displaystyle\tan\theta(k)\Bigr{)}\tau_{c}(k)\biggl{.}-r^{2}\cos^{2}\theta(k)\Bigl{(}\frac{\alpha_{1}(k)}{l_{1}(k)}+\frac{\alpha_{2}(k)}{l_{2}(k)}\Bigr{)}\biggr{]},$
(4a) $\displaystyle\bar{F}_{2}(k)=$ $\displaystyle\
\bar{F}_{1}(k)-\frac{\tau_{c}(k)}{r\cos\theta(k)}.$ (4b)

####  II-C2 Pressure into Contraction Force

The pressure-to-contraction force block shown in Fig. 2 estimates the
contraction force of the PAM by using the following experimental polynomial
function of the inner pressure and measured joint angle:

$\displaystyle\hat{F}_{i}$
$\displaystyle(\theta(k),P_{i}(k))=\hat{a}_{i}(\theta(k))P_{i}(k)+\hat{b}_{i}(\theta(k)),$
(5a)
$\displaystyle\hat{a}_{i}(\theta)=p^{\hat{a}_{1}}_{i}\theta+p^{\hat{a}_{2}}_{i},\quad\hat{b}_{i}(\theta)=p^{\hat{b}_{1}}_{i}\theta+p^{\hat{b}_{2}}_{i},$
(5b)

where $\hat{F}_{i}$, $\hat{a}_{i}$, and $\hat{b}_{i}$, $\forall
i\in\mathcal{I}$, are an estimate of the PAM contraction force and the
coefficients of (5a), respectively.

A reason for introducing the contraction force model (5) in addition to (1) is
to reduce the number of approximation, which will facilitate the discussion in
Section III. Model (1) involves a nonlinear term $\sin\theta$, which requires
the approximation process to be a polynomial. On the other hand, the model (5)
is a polynomial of $\theta$ and captures a static relationship between the
angle and the contraction force. Indeed, the fitting results of (5) are shown
in Fig. 3. In Fig. 3LABEL:sub@sfig:P-F, the black circles represent
experimentally measured data, and the colored lines represent fitting results
using (5) for each angle. In this case, the coefficients (5b) for each angle
correspond to Fig. 3LABEL:sub@sfig:v-w, where the black circles represent the
coefficients corresponding to the colored lines in Fig. 3LABEL:sub@sfig:P-F.
The fitting test confirms that the $\theta$-dependent model (5) is valid.

(a)

(b)

Figure 3: Fitting result of the $\theta$-dependent model (5) for capturing
static characteristics with regard to $P$, $F$, and $\theta$ of the PAMs. (a)
$P$-$F$ relationship of the PAM. (b) $\hat{a}$-$\theta$ and $\hat{b}$-$\theta$
relationships of the PAM.

####  II-C3 PI Controller

The proposed controller includes two types of PI controllers to compensate for
errors in the joint angle and the contraction force. The controller for the
joint angle uses a feedback error between $\bar{\theta}$ and $\theta$, denoted
by $e^{\theta}:=\bar{\theta}-\theta$, to compute the command torque
$\tau_{c}$, as follows:

$\displaystyle\begin{bmatrix}x^{\theta}(k+1)\\\
\tau_{c}(k)\end{bmatrix}=\begin{bmatrix}1&T_{s}\\\ G_{\rm I}^{\theta}&G_{\rm
P}^{\theta}\end{bmatrix}\begin{bmatrix}x^{\theta}(k)\\\
e^{\theta}(k)\end{bmatrix},$ (6)

where $x^{\theta}$ is the state and $G_{\rm P}^{\theta}$ and $G_{\rm
I}^{\theta}$ are the proportional and integral gains, respectively.

Regarding the contraction force controller, the outputs of the reference
generator, $\bar{F}_{1}$ and $\bar{F}_{2}$, are compared with the estimated
forces $\hat{F_{1}}$ and $\hat{F_{2}}$ computed by (5), respectively. The
resulting errors are fed to the PI controllers to generate the control input
voltages $u_{1}$ and $u_{2}$ for each PDCV. The controllers are given as
follows:

$\displaystyle\begin{bmatrix}x_{i}^{F}(k+1)\\\
u_{i}(k)\end{bmatrix}=\begin{bmatrix}1&T_{s}&0\\\ G_{\rm I}^{F}&G_{\rm
P}^{F}&1\end{bmatrix}\begin{bmatrix}x_{i}^{F}(k)\\\ e_{i}^{F}(k)\\\
\beta_{i}\end{bmatrix},\ \forall i\in\mathcal{I},$ (7)

where $x_{i}^{F},\forall i\in\mathcal{I},$ is the state of the $i$-th
controller, $e_{i}^{F}:=\bar{F}_{i}-\hat{F}_{i},\forall i\in\mathcal{I},$ is
the error in the contraction force of the $i$th PAM, $G_{\rm P}^{F}$ and
$G_{\rm I}^{F}$ are the proportional and integral gains, respectively, and the
gains are common between the controllers. Additionally, we set
$\beta_{i}=5.0,\forall i\in\mathcal{I}$.

Consequently, defining the controller state by $x:=[\,x^{\theta}\ \,x_{1}^{F}\
\,x_{2}^{F}\,]^{\mathrm{T}}$ and eliminating $\tau_{c}$, $e^{\theta}$,
$e^{F}_{i}$, $\bar{F}_{i}$, and $\hat{F}_{i}$, $\forall i\in\mathcal{I}$ from
(4)–(7) results in the nonlinear controller (3) with the coefficients and
functions in (8).

$\displaystyle A(k;\theta)=$ $\displaystyle\begin{bmatrix}1&0&0\\\
\frac{T_{s}G_{\textrm{I}}^{\theta}l_{1}(k)l_{2}(k)}{r^{2}(l_{1}(k)+l_{2}(k))\cos^{2}\theta(k)}\Bigl{(}\frac{r\cos\theta(k)}{l_{2}(k)}-\tan\theta(k)\Bigr{)}&1&0\\\
T_{s}G_{\textrm{I}}^{\theta}\left\\{\frac{l_{1}(k)l_{2}(k)}{r^{2}(l_{1}(k)+l_{2}(k))\cos^{2}\theta(k)}\Bigl{(}\frac{r\cos\theta(k)}{l_{2}(k)}-\tan\theta(k)\Bigr{)}-\frac{1}{r\cos\theta(k)}\right\\}&0&1\end{bmatrix},\quad
g(k;\zeta)=T_{s}\begin{bmatrix}\bar{\theta}(k)-\theta(k)\\\ h_{1}(k;\zeta)\\\
h_{2}(k;\zeta)\end{bmatrix},$ (8a) $\displaystyle C(k;\theta)=$
$\displaystyle\begin{bmatrix}G_{\textrm{P}}^{F}G_{\textrm{I}}^{\theta}&G_{\textrm{I}}^{F}&0\\\
G_{\textrm{P}}^{F}G_{\textrm{I}}^{\theta}\left\\{\frac{l_{1}(k)l_{2}(k)}{r^{2}(l_{1}(k)+l_{2}(k))\cos^{2}\theta(k)}\Bigl{(}\frac{r\cos\theta(k)}{l_{2}(k)}-\tan\theta(k)\Bigr{)}-\frac{1}{r\cos\theta(k)}\right\\}&0&G_{I}^{F}\end{bmatrix},\quad
h(k;\zeta)=G_{\textrm{P}}^{F}\begin{bmatrix}h_{1}(k;\zeta)\\\
h_{2}(k;\zeta)\end{bmatrix},$ (8b) $\displaystyle h_{1}(k;\zeta):=$
$\displaystyle\frac{l_{1}(k)l_{2}(k)}{l_{1}(k)+l_{2}(k)}\left[\frac{1}{r^{2}\cos^{2}\theta(k)}\left\\{-\bar{K}_{P}(k)+G_{\textrm{P}}^{\theta}\Bigl{(}\frac{r\cos\theta(t)}{l_{2}(k)}-\tan\theta(k)\Bigr{)}(\bar{\theta}(k)-\theta(k))\right\\}+\frac{\alpha_{1}(k)}{l_{1}(k)}+\frac{\alpha_{2}(k)}{l_{2}(k)}\right]-$
$\displaystyle\hat{a}_{1}P_{1}(k)-\hat{b}_{1},$ (8c) $\displaystyle
h_{2}(k;\zeta):=$
$\displaystyle\frac{l_{1}(k)l_{2}(k)}{l_{1}(k)+l_{2}(k)}\left[\frac{1}{r^{2}\cos^{2}\theta(k)}\left\\{-\bar{K}_{P}(k)+G_{\textrm{P}}^{\theta}\Bigl{(}\frac{r\cos\theta(k)}{l_{2}(k)}-\tan\theta(k)-\frac{l_{1}(k)+l_{2}(k)}{l_{1}(k)l_{2}(k)}r\cos\theta(k)\Bigr{)}(\bar{\theta}(k)-\theta(k))\right\\}+\right.$
$\displaystyle\left.\frac{\alpha_{1}(k)}{l_{1}(k)}+\frac{\alpha_{2}(k)}{l_{2}(k)}\right]-\hat{a}_{2}P_{2}(k)-\hat{b}_{2},$
(8d)

(a)

(b)

(c)

(d)

Figure 4: Simulation results of the original simultaneous control of joint
angle and stiffness for Reference 2. (a) Joint angle and stiffness, (b)
tracking errors of the joint angle and stiffness, (c) control inputs to the
valve, and (d) inner pressure of each PAM.

###  II-D Numerical Verification

This section presents the numerical simulations conducted to evaluate the
proposed nonlinear controller in terms of the control performance with respect
to simultaneous angle-stiffness control. For the simulations, we used the
following state-space model of an antagonistic PAM system [40]:
$x_{p}(k+1)=f_{\sigma}(x_{p}(k),u(k))$ if $x_{p}(k)\in\mathcal{X}_{\sigma}$
holds, $y(k)=h_{p}(x_{p}(k))$, where $u:=[u_{1}\
u_{2}]^{\mathrm{T}}\in[0,10]^{2}\subset\mathbb{R}^{2}$ is the input voltages
(V) to the two PDCVs, the state variable is $x_{p}:=[\,\theta\,\
\dot{\theta}\,\ P_{1}\,\ P_{2}\,]^{\mathrm{T}}\in\mathbb{R}^{4}$, and the
output variable is $y:=[\,\theta\,\ P_{1}\,\ P_{2}\,\
K_{P}\,]^{\mathrm{T}}\in\mathbb{R}^{4}$, in which a set of allowable absolute
pressures (kPa) is determined by the specification of the PAMs, i.e.,
$P_{1},P_{2}\in[200,750]$,
$f_{\sigma}:\mathbb{R}^{4}\rightarrow\mathbb{R}^{4}$ is a nonlinear function
with 18 subsystems, and it switches according to if-then rules,
${\mathcal{X}}_{\sigma}:=\\{x_{p}\in\mathbb{R}^{4}\,|\,\Psi_{\sigma}(x_{p})>0\\}$
is the state set, where $\sigma\in\\{1,2,\cdots,18\\}$ is the index of the
subsystem, $\Psi_{\sigma}(x_{p})$ is a function derived from the modes in the
form of if-then rules, and the function
$h_{p}:\mathbb{R}^{4}\rightarrow\mathbb{R}^{4}$ is an observation equation.

We considered two types of references of the joint angle and stiffness.
Reference 1 was set to

$\displaystyle(\bar{\theta}(k),\bar{K}_{P}(k))=$ $\displaystyle\
\begin{cases}\,(10,8)&\textrm{if}\ \ 0\leq T_{s}k<15,\\\ \,(10,6)&\textrm{if}\
\ 15\leq T_{s}k<30,\\\ \,(10,4)&\textrm{if}\ \ 30\leq T_{s}k<45,\\\
\end{cases}$

and Reference 2 was set to

$\displaystyle(\bar{\theta}(k),\bar{K}_{P}(k))=$ $\displaystyle\
\begin{cases}\,(5,9)&\textrm{if}\ \ 0\leq T_{s}k<15,\\\ \,(15,6)&\textrm{if}\
\ 15\leq T_{s}k<30,\\\ \,(10,7)&\textrm{if}\ \ 30\leq T_{s}k<45.\\\
\end{cases}$

The parameters of the PI controllers in (6) and (7) were determined by trial
and error, resulting in $G_{\rm P}^{\theta}=0.25$, $G_{\rm I}^{\theta}=0.13$,
$G_{\rm P}^{F}=0.088$, and $G_{\rm I}^{F}=0.08$. The initial voltage commands
$u_{1}(0)=u_{2}(0)=5.5$ were given before the start of the control, and the
control was started after sufficient time had elapsed. Additionally, note that
the simulation and experimental results regarding Reference 1 have been
omitted, and the results are discussed in Section IV-C owing to page
limitations.

The simulation results are provided in Fig. 4, where (a) shows the time
response of the joint angle and stiffness, (b) presents the tracking errors of
the joint angle and stiffness, (c) depicts the control inputs to the valve,
and (d) shows the inner pressure of each PAM. Figs.
4LABEL:sub@sfig:ori2_outLABEL:sub@sfig:ori2_error confirm that the joint angle
and stiffness track the reference in a steady state. In this case, the input
voltage oscillates in a transient state, as shown in Fig.
4LABEL:sub@sfig:ori2_ele. The operating pressure behavior is shown in Fig.
4LABEL:sub@sfig:ori2_pre. These results confirm that the proposed nonlinear
controller achieves simultaneous tracking control of the joint angle and
stiffness of the PAM actuator system.

However, the coefficients and input functions in the proposed controller are
rational functions, which makes them difficult to implement in a controller
device in an encrypted control fashion. Therefore, in the next section, we
consider polynomial approximation of the proposed controller.

###### Remark 1

The novelty of the proposed controller (3) is the use of the
$\theta$-dependent contraction force model (5), compared to the UKF-based
simultaneous controller [24]. The UKF-based controller involving fluid
dynamics [24] and sliding mode controllers [20, 21] include an if-then rule;
therefore, obtaining an alternative controller in a polynomial is difficult.
Meanwhile, the nonlinear controllers proposed in [19, 22, 23] can be applied
to a polynomial approximation, whereas costly force sensors are required to
measure the contraction force.

##  III Polynomial Approximation of The Controller

This section presents the approximation of the proposed controller into
polynomial functions and investigates the impact of the approximation on the
response of the experimental PAM control system.

###  III-A Controller Approximation

This study uses one- and two-variable functions to approximate the reference
generator (4) in polynomials because only the reference generator is a
rational function. By introducing five functions denoted by
$f_{1},\dots,f_{5}$, let the rational functions (4) be rewritten as follows:

$\displaystyle\bar{F}_{1}(k)=$
$\displaystyle\,f_{1}(k;\theta,\bar{K}_{P})+f_{2}(k;\theta)\tau_{c}(k)+f_{3}(k;\theta,P_{1})+$
$\displaystyle f_{4}(k;\theta,P_{2}),$ $\displaystyle\bar{F}_{2}(k)=$
$\displaystyle\,\bar{F}_{1}(k)+f_{5}(k;\theta)\tau_{c}(k),$

where

$\displaystyle f_{1}(k;\theta,\bar{K}_{P}):=$
$\displaystyle\,-\frac{l_{1}(k)l_{2}(k)\bar{K}_{P}(k)}{r^{2}(l_{1}(k)+l_{2}(k))\cos^{2}\theta(k)},$
$\displaystyle f_{2}(k;\theta):=$
$\displaystyle\,\frac{l_{1}(k)l_{2}(k)}{r^{2}(l_{1}(k)+l_{2}(k))\cos^{2}\theta(k)}\Bigl{(}\frac{r\cos\theta(k)}{l_{2}(k)}-$
$\displaystyle\,\tan\theta(k)\Bigr{)},$ $\displaystyle
f_{3}(k;\theta,P_{1}):=$
$\displaystyle\,\frac{\alpha_{1}(P_{1}(k))l_{2}(k)}{l_{1}(k)+l_{2}(k)},$
$\displaystyle f_{4}(k;\theta,P_{2}):=$
$\displaystyle\,\frac{\alpha_{2}(P_{2}(k))l_{1}(k)}{l_{1}(k)+l_{2}(k)},$
$\displaystyle f_{5}(k;\theta):=$ $\displaystyle\,-\frac{1}{r\cos\theta(k)}.$

For the five functions, we used the LASSO-based polynomial approximation [41]
and manually removed coefficients with relatively small values, which resulted
in the following approximated functions:

$\displaystyle f_{1}(k;\theta,\bar{K}_{P})\approx$
$\displaystyle\,w_{1}\bar{K}_{p}(k)+w_{2}\theta^{2}(k)\bar{K}_{P}(k)+w_{3}\,=:\hat{f_{1}},$
$\displaystyle f_{2}(k;\theta)\approx$
$\displaystyle\,w_{4}\theta(k)+w_{5}\theta^{2}(k)+w_{6}\,=:\hat{f_{2}},$
$\displaystyle f_{3}(k;\theta,P_{1})\approx$
$\displaystyle\,w_{7}P_{1}(k)+w_{8}\theta(k)P_{1}(k)+$
$\displaystyle\,w_{s}\theta(k)P_{1}^{2}(k)+w_{9},$ $\displaystyle\approx$
$\displaystyle\,w_{7}P_{1}(k)+w_{8}\theta(k)P_{1}(k)+w_{9}\,=:\hat{f_{3}},$
$\displaystyle f_{4}(k;\theta,P_{2})\approx$
$\displaystyle\,w_{10}P_{2}(k)+w_{11}\theta(k)P_{2}(k)+w_{12}\,=:\hat{f_{4}},$
$\displaystyle f_{5}(k;\theta)\approx$
$\displaystyle\,w_{13}\theta^{2}(k)+w_{14}\,=:\hat{f_{5}},$

where $w_{i},\,\forall i\in\\{1,2,\cdots,14\\}$ and $w_{s}$ are coefficients
of approximated polynomials $\hat{f}_{j},\,\forall j\in\\{1,2,\cdots,5\\}$,
and their values are listed in TABLE II, where the regularization parameter of
LASSO was set to 1.0. $\hat{f_{2}}$ and $\hat{f_{5}}$ are cubic polynomials in
$\theta$, and the other functions are cubic in both arguments. Furthermore,
$w_{s}$ of $\hat{f}_{3}$ resulted in a nonzero term, but it was significantly
small compared to the other coefficients; thus we removed $w_{s}$.

Consequently, we obtained the polynomial controller in (3) by using the
approximated coefficient matrices and functions specified below,

$\displaystyle A(k;\theta)\approx$ $\displaystyle\,\begin{bmatrix}1&0&0\\\
T_{s}G_{\textrm{I}}^{\theta}\hat{f_{2}}(k;\theta)&1&0\\\
T_{s}G_{\textrm{I}}^{\theta}(\hat{f_{2}}(k;\theta)+\hat{f_{5}}(k;\theta))&0&1\end{bmatrix},$
(9a) $\displaystyle C(k;\theta)\approx$
$\displaystyle\,\begin{bmatrix}G_{\textrm{P}}^{F}G_{\textrm{I}}^{\theta}\hat{f_{2}}(k;\theta)&G_{\textrm{I}}^{F}&0\\\
G_{\textrm{P}}^{F}G_{\textrm{I}}^{\theta}(\hat{f_{2}}(k;\theta)+\hat{f_{5}}(k;\theta))&0&G_{\textrm{I}}^{F}\end{bmatrix},$
(9b) $\displaystyle h_{1}(k;\zeta)\approx$
$\displaystyle\,\hat{f_{1}}(k;\theta,\bar{K}_{P})+G_{\textrm{P}}^{\theta}\hat{f_{2}}(k;\theta)(\bar{\theta}(k)-\theta(k))+$
$\displaystyle\,\hat{f_{3}}(k;\theta,P_{1})+\hat{f_{4}}(k;\theta,P_{2})-\hat{a}_{1}P_{1}(k)-\hat{b}_{1},$
(9c) $\displaystyle h_{2}(k;\zeta)\approx$
$\displaystyle\,\hat{f_{1}}(k;\theta,\bar{K}_{P})+\hat{f_{3}}(k;\theta,P_{1})+\hat{f_{4}}(k;\theta,P_{2})+$
$\displaystyle\,G_{\textrm{P}}^{\theta}(\hat{f_{2}}(k;\theta)+\hat{f_{5}}(k;\theta))(\bar{\theta}(k)-\theta(k))-$
$\displaystyle\,\hat{a}_{2}P_{2}(k)-\hat{b}_{2}.$ (9d)

TABLE II: List of Parameters in Approximated Controller $w_{1}$ | $-61.7$ | $p_{1}^{a_{1}}$ | 7.05$\times 10^{-3}$
---|---|---|---
$w_{2}$ | $-1.89\times 10^{-2}$ | $p_{1}^{a_{2}}$ | $-1.02\times 10^{-4}$
$w_{3}$ | $-1.58$ | $p_{1}^{b_{1}}$ | $-5.57\times 10^{2}$
$w_{4}$ | $13.6$ | $p_{1}^{b_{2}}$ | $72.86$
$w_{5}$ | $-1.24$ | $p_{2}^{a_{1}}$ | $6.42\times 10^{3}$
$w_{6}$ | $1.73\times 10^{-3}$ | $p_{2}^{a_{2}}$ | $-9.18\times 10^{-4}$
$w_{7}$ | $-5.12\times 10^{-1}$ | $p_{2}^{b_{1}}$ | $-1.98\times 10^{2}$
$w_{8}$ | $-1.18\times 10^{-3}$ | $p_{2}^{b_{2}}$ | 15.75
$w_{s}$ | $-5.76\times 10^{-7}$ | $p_{1}^{\hat{a}_{1}}$ | $-2.55\times 10^{-4}$
$w_{9}$ | 36.5 | $p_{1}^{\hat{a}_{2}}$ | $1.45\times 10^{-4}$
$w_{10}$ | $-4.31\times 10^{-1}$ | $p_{1}^{\hat{b}_{1}}$ | 20.14
$w_{11}$ | $1.54\times 10^{-3}$ | $p_{1}^{\hat{b}_{2}}$ | $-19.01$
$w_{12}$ | $-9.35$ | $p_{2}^{\hat{a}_{1}}$ | $2.18\times 10^{-4}$
$w_{13}$ | $-27.3$ | $p_{2}^{\hat{a}_{2}}$ | $1.45\times 10^{-4}$
$w_{14}$ | $-3.92\times 10^{-3}$ | $p_{2}^{\hat{b}_{1}}$ | $1.35$
$G_{\textrm{P}}^{\theta}$ | 1.3 | $p_{2}^{\hat{b}_{2}}$ | $-18.67$
$G_{\textrm{I}}^{\theta}$ | $2.43\times 10^{-1}$ | $G_{\textrm{I}}^{F}$ | $2.50\times 10^{-2}$
$G_{\textrm{P}}^{F}$ | $8.80\times 10^{-2}$ |  |

###  III-B Evaluation of Approximation Impacts

The impact of the controller approximation on the control performance is
evaluated by comparing the experimental time responses of the control systems
using the original controller (8) and the approximated controller (9). The
scenarios of the control experiments are simultaneous tracking control of the
joint angle and stiffness with and without a load. The references of the joint
angle and stiffness set in the experiments are the same as those employed in
the simulations. Before the start of the control, an initial voltage command
$u_{1}(0)=u_{2}(0)=5.5$ were provided, and the control was initiated after
sufficient time had elapsed.

The experimental results for the original and approximated controllers without
a load are shown in Fig. 5. In this figure, the blue and red lines indicate
the results of the original and approximated control systems, respectively,
and their meanings are the same as those in Fig. 4. Figs.
5LABEL:sub@sfig:app2_outLABEL:sub@sfig:app2_error confirm that the angle and
stiffness track the reference in the steady state and that sufficient control
performance is maintained. Additionally, the improvement in transient-response
control performance is believed to be due to the impacts of quantization;
however, the details have not been clearly established. Figs.
5LABEL:sub@sfig:app2_eleLABEL:sub@sfig:app2_pre observe that the experimental
results of the original and approximated controllers exhibit similar behavior
in the steady state.

Moreover, we conducted a control experiment in which a 1.5 kg load was hung
from the left side of the joint before starting control to evaluate the impact
of the approximation. The experimental results of the loads are shown in Fig.
6, where the line colors in each figure are the same as those in Fig. 5. Figs.
6LABEL:sub@sfig:app4_outLABEL:sub@sfig:app4_error confirm that the angle and
stiffness track the reference in the steady state. Good control tracking
causes the controller to compensate for the impact of the load, which can be
observed as a difference of approximately 50 kPa over the steps between the
responses of the inner pressures, as shown in Figs. 5LABEL:sub@sfig:app2_pre
and 6LABEL:sub@sfig:app4_pre. Similarly, Fig. 6LABEL:sub@sfig:app4_ele
demonstrates that the experimental results of the two controllers exhibit
similar behaviors in the steady state, respectively.

The experimental results confirm that the approximated controller achieves
almost the same control performance as the original controller, implying that
the impact of the controller approximation is negligible. Therefore, the
original controller can be replaced with an approximated controller that is
adequate for secure implementation.

(a)

(b)

(c)

(d)

Figure 5: Comparison of experimental results of original and approximated
polynomial load-free controls of joint angle and stiffness for Reference 2.
(a) Joint angle and stiffness, (b) tracking errors of the joint angle and
stiffness, (c) control inputs to the valve, and (d) inner pressure of each
PAM.

(a)

(b)

(c)

(d)

Figure 6: Comparison of experimental results of original and approximated
polynomial controls with load for Reference 2. (a) Joint angle and stiffness,
(b) Tracking errors of the joint angle and stiffness, (c) Control inputs to
the valve, and (d) Inner pressure of each PAM.

##  IV Secure PAM Actuator with Encrypted Polynomial-Type Controller

This section describes the encrypted control of the joint angle and the
stiffness of an antagonistic PAM actuator system.

###  IV-A Secure Implementation

The controller encryption method in [33] can be applied to a linear controller
in a matrix-vector product form. Hence, we consider to transform our proposed
polynomial controller into a product of matrix and vector. In this study, by
defining a vector of monomials as $\xi\in\mathbb{R}^{18}$, the polynomial
controller (3) with the coefficients and functions (9) can be transformed into
the following matrix-vector multiplication:

$\displaystyle\psi(k)=$ $\displaystyle\,\Phi\xi(k),\ \,\forall
k\in\mathbb{Z}^{+},$ (10)

with

$\displaystyle\xi:=$ $\displaystyle\,[\,\bar{K}_{P}\,\
\bar{K}_{P}\theta^{2}\,\ \bar{\theta}\,\ \theta\,\ x^{\theta}\,\
\bar{\theta}\theta\,\ \theta^{2}\,\ x^{\theta}\theta\,\
\bar{\theta}\theta^{2}\,\ \theta^{3}\,\ x^{\theta}\theta^{2}\,\ P_{1}$
$\displaystyle\ \ \theta P_{1}\,\ P_{2}\,\ \theta P_{2}\,\ 1\,\ x^{F_{1}}\,\
x^{F_{2}}\,]^{\mathrm{T}},$ $\displaystyle\psi(k):=$
$\displaystyle\,[\,x^{\theta}(k+1)\,\ x^{F}_{1}(k+1)\,\ x^{F}_{2}(k+1)\,\
u_{1}(k)\,\ u_{2}(k)\,]^{\mathrm{T}},$

where $\xi$ and $\psi\in\mathbb{R}^{5}$ are the input and output of the
encrypted controller, respectively. The coefficient matrix
$\Phi\in\mathbb{R}^{5\times 18}$ is given as follows:

$\displaystyle\Phi=\left[\begin{array}[]{cccccc}0&0&T_{s}&-T_{s}&1&0\\\
T_{s}w_{1}&T_{s}w_{2}&T_{s}\phi_{3}&T_{s}\phi_{4}&T_{s}\phi_{5}&T_{s}\phi_{6}\\\
T_{s}w_{1}&T_{s}w_{2}&T_{s}\nu_{3}&T_{s}\nu_{4}&T_{s}\nu_{5}&T_{s}\nu_{6}\\\
G_{P}^{F}w_{1}&G_{P}^{F}w_{2}&G_{P}^{F}\phi_{3}&G_{P}^{F}\phi_{4}&G_{P}^{F}\phi_{5}&G_{P}^{F}\phi_{6}\\\
G_{P}^{F}w_{1}&G_{P}^{F}w_{2}&G_{P}^{F}\nu_{3}&G_{P}^{F}\nu_{4}&G_{P}^{F}\nu_{5}&G_{P}^{F}\nu_{6}\\\
\end{array}\right.$ $\displaystyle\left.\begin{array}[]{cccccc}0&0&0&0&0&0\\\
T_{s}\phi_{7}&T_{s}\phi_{8}&T_{s}\phi_{9}&T_{s}\phi_{10}&T_{s}\phi_{11}&T_{s}\phi_{12}\\\
G_{P}^{F}\phi_{7}&G_{P}^{F}\phi_{8}&G_{P}^{F}\phi_{9}&G_{P}^{F}\phi_{10}&G_{P}^{F}\phi_{11}&T_{s}\nu_{12}\\\
G_{P}^{F}\phi_{7}&G_{P}^{F}\phi_{8}&G_{P}^{F}\phi_{9}&G_{P}^{F}\phi_{10}&G_{P}^{F}\phi_{11}&G_{P}^{F}\phi_{12}\\\
G_{P}^{F}\nu_{7}&G_{P}^{F}\nu_{8}&G_{P}^{F}\nu_{9}&G_{P}^{F}\nu_{10}&G_{P}^{F}\nu_{11}&G_{P}^{F}\nu_{12}\\\
\end{array}\right.$ $\displaystyle\left.\begin{array}[]{cccccc}0&0&0&0&0&0\\\
T_{s}\phi_{13}&T_{s}\phi_{14}&T_{s}\phi_{15}&T_{s}\phi_{16}&1&0\\\
T_{s}\nu_{13}&T_{s}\nu_{14}&T_{s}\nu_{15}&T_{s}\nu_{16}&0&1\\\
G_{P}^{F}\phi_{13}&G_{P}^{F}\phi_{14}&G_{P}^{F}\phi_{15}&G_{P}^{F}\phi_{16}+\beta_{1}&G_{I}^{F}&0\\\
G_{P}^{F}\nu_{13}&G_{P}^{F}\nu_{14}&G_{P}^{F}\nu_{15}&G_{P}^{F}\nu_{16}+\beta_{2}&0&G_{I}^{F}\end{array}\right],$

where $\phi_{3}:=w_{4}G_{P}^{\theta}$,
$\phi_{4}:=-p_{1}^{b_{1}}-w_{4}G_{P}^{\theta}$,
$\phi_{5}:=w_{4}G_{I}^{\theta}$, $\phi_{6}:=w_{5}G_{P}^{\theta}$,
$\phi_{7}:=-w_{5}G_{P}^{\theta}$, $\phi_{8}:=w_{5}G_{I}^{\theta}$,
$\phi_{9}:=w_{6}G_{P}^{\theta}$, $\phi_{10}:=-w_{6}G_{P}^{\theta}$,
$\phi_{11}:=w_{6}G_{I}^{\theta}$, $\phi_{12}:=w_{7}-p_{2}^{a_{1}}$,
$\phi_{13}:=w_{8}-p_{1}^{\hat{a}_{1}}$, $\phi_{14}:=w_{10}$,
$\phi_{15}:=w_{11}$, $\phi_{16}:=w_{3}+w_{9}+w_{12}-p_{2}^{\hat{b}_{1}}$,
$\nu_{3}:=(w_{4}+w_{13})G_{P}^{\theta}$,
$\nu_{4}:=-p_{1}^{\hat{b}_{2}}-(w_{4}+w_{13})G_{P}^{\theta}$,
$\nu_{5}:=(w_{4}+w_{13})G_{I}^{\theta}$, $\nu_{6}:=w_{5}G_{P}^{\theta}$,
$\nu_{7}:=-w_{5}G_{P}^{\theta}$, $\nu_{8}:=w_{5}G_{I}^{\theta}$,
$\nu_{9}:=(w_{6}+w_{14})G_{P}^{\theta}$,
$\nu_{10}:=-(w_{6}+w_{14})G_{P}^{\theta}$,
$\nu_{11}:=(w_{6}+w_{14})G_{I}^{\theta}$, $\nu_{12}:=w_{7}$,
$\nu_{13}:=w_{8}$, $\nu_{14}:=w_{10}-p_{2}^{\hat{a}_{2}}$,
$\nu_{15}:=w_{11}-p_{1}^{\hat{a}_{2}}$, and
$\nu_{16}:=w_{3}+w_{9}+w_{12}-p_{2}^{\hat{b}_{2}}$. The augmented formulation
of the controller (10) is a linear operation; thus, the same procedure for
encrypting a linear controller is applied to encrypt the controller, as
described in [36]. The configuration of the control system is illustrated in
Fig. 7, where the controller input $\xi$ is generated by the reference
signals, measurements, and states of the PI controllers. The input $\xi$ is
encrypted by using an ElGamal encryption scheme in the Enc block. The
encrypted controller $\mathsf{Enc}(\Phi)$ conducts multiplicative homomorphic
operations, and $\textsf{Dec}^{+}$ extracts the plaintext states and control
inputs to each PDCV. For the notations regarding the encryption scheme, please
refer to [33].

Figure 7: Encrypted angle and stiffness control system.

###  IV-B Experimental Validation

This section verifies the encrypted control system by comparing the
experimental results obtained using the encrypted and original control
systems. A key length of 64 bit was chosen, and the scaling parameters
$\Delta_{\xi}$ and $\Delta_{\Phi}$ were set to $1.0\times 10^{8}$, which were
introduced in [33] to manage the quantization errors due to encryption.

The scenarios of the control experiments involve simultaneous tracking control
of the joint angle and stiffness with and without the load. The references of
the joint angle and stiffness were the same as those of the simulations in
Figs. 4. Before the start of the control, an initial voltage command
$u_{1}(0)=u_{2}(0)=5.5$ was provided, and the control was initiated after
sufficient time had elapsed. The experimental results are shown in Figs. 10.
In the figures, the blue and green lines represnet the results of the original
and encrypted control systems, respectively. Figs.
9LABEL:sub@sfig:enc2_outLABEL:sub@sfig:enc2_error and
10LABEL:sub@sfig:enc4_outLABEL:sub@sfig:enc4_error confirm that the angle and
stiffness track the reference in the steady state and that sufficient control
performance is maintained. Similarly,they show that the controller compensated
for the impact of the load, which can be observed as a difference of
approximately 50 kPa.

The computation time involved in the encrypted control system is shown in Fig.
8, where Figs. 8LABEL:sub@fig:cal_time_bLABEL:sub@fig:cal_time_d correspond to
Figs. 9 and 10. The figures show that the control operation is in real time
because the computation time for each step is less than the sampling period of
20 ms.

(a)

(b)

Figure 8: Processing time: less than the sampling period of 20 ms. (a)
Reference 2 without load and (b) Reference 2 with load.

(a)

(b)

(c)

(d)

Figure 9: Experimental results of unencrypted and encrypted simultaneous load-
free control of joint angle and stiffness for Reference 2. (a) Joint angle and
stiffness, (b) tracking errors of the joint angle and stiffness, (c) control
inputs to the valve, and (d) inner pressure of each PAM.

(a)

(b)

(c)

(d)

Figure 10: Experimental results of unencrypted and encrypted simultaneous
controls with load for Reference 2. (a) Joint angle and stiffness, (b)
tracking errors of the joint angle and stiffness, (c) control inputs to the
valve, and (d) inner pressure of each PAM.

###  IV-C Quantitative Investigation

We quantitatively evaluate the control results presented so far to investigate
the impact of the secure implementation of the nonlinear controller on the
control performance. For this purpose, we introduce the $\ell_{2}$ norm to
measure the tracking error during a specific duration for each control result,
as shown in Figs. 5-6 and 9-10. The $\ell_{2}$ norm is defined as
$\gamma(z,\bar{z}):=\sqrt{\sum_{k=k_{0}}^{k_{1}}(z(k)-\bar{z}(k))^{2}}$, where
$z$ and $\bar{z}$ represent scalar signal sequences, and the evaluation steps
(closed interval) $[k_{0},k_{1}]$ are set to $[500,749]$, $[1250,1499]$, and
$[2000,2249]$. These intervals are labeled #1, #2, and #3, respectively, in
References 1 and 2.

The $\ell_{2}$-norm scores of the original, approximated, and encrypted
controls with and without loads are summarized in Figs. 11 and 12,
respectively. In each figure, LABEL:sub@sfig:con1_angle and
LABEL:sub@sfig:con1_stiff display $\gamma(\theta,\bar{\theta})$ and
$\gamma(K_{P},\bar{K}_{P})$, respectively, for the three evaluation steps. The
blue, red, and green colors represent the $\ell_{2}$-norm values of the
original, approximated, and encrypted controllers, respectively. The dots and
error bars indicate the average, maximum, and minimum $\ell_{2}$-norm values
from 10 experiments.

These figures confirm that the resulting scores tend to increase as the
procedure advances through the polynomial approximation and secure
implementation. However, in many cases, the scores of the original and
encrypted controllers are similar. The worst rate of the controlled signal
relative to its reference is $3.89/4.00\times$ 100=97.3 % for #3 of Reference
1 in Fig. 11LABEL:sub@sfig:con1_stiff. This result implies that the proposed
encrypted control achieved a tracking error of less than 2.7 %. Moreover, in
cases where a relatively large change in the score occurs, such as #2 of
Reference 2 in Fig. 11LABEL:sub@sfig:con1_angle and #3 of Reference 1 in Fig.
12LABEL:sub@sfig:con2_stiff_load, the effect of the polynomial approximation
prevails over the change in. In other words, the change in score from blue to
red is larger than that from red to green. This discussion provides further
insight, suggesting that increasing the accuracy of the polynomial
approximation could help avoid the degradation of the original control
performance. It is crucial to consider an accuracy-aware approximation method
suitable for secure implementation, which will be the focus of future work. In
addition, these figures confirm that, owing to the lack of significant
differences between the cases without and with the load, the proposed control
system is capable of compensating for the unknown load.

(a)

(b)

Figure 11: $\ell_{2}$-norm scores of the tracking errors of the (a) joint
angle and (b) stiffness using the original, approximated, and encrypted
controllers without a load.

(a)

(b)

Figure 12: $\ell_{2}$-norm scores of the tracking error of the (a) joint angle
and (b) stiffness using the original, approximated, and encrypted controllers
with a load.

##  V Conclusion

This study proposed an encrypted simultaneous control system for antagonistic
PAM actuators aimed at developing cyber-secure and safe PAM actuator systems.
A novel nonlinear controller was designed to track the joint angle and
stiffness simultaneously based on the PAM actuator model. By applying the
polynomial approximation technique to a nonlinear controller, we obtained a
polynomial-type controller. Subsequently, through the secure implementation in
the control device, we developed a secure PAM actuator system. For
experimental validation, the $\ell_{2}$ norm was introduced to measure and
compare the experimental results of the original, approximated, and encrypted
controllers. The experimental results showed that the proposed encrypted
controller achieved simultaneous tracking of the joint angle and stiffness of
the PAM actuator with a tracking error of less than 2.7 %. Consequently, the
developed PAM actuator system, enabled by the secure implementation of the
simultaneous controller, enhances security while maintaining a control
performance similar to that of the original controller. The developed actuator
is expected to be used in secure and safe PAM-driven devices, such as nursing
care robots, rehabilitation orthoses, and power-assisted orthoses for remote
usage applications.

In future work, we plan to improve the control performance of the encrypted
controller further by addressing the nonlinear characteristics of PAMs,
including Coulomb friction and fluid dynamics. To mitigate the performance
degradation caused by the polynomial approximation, we investigate more
effective strategies for tuning the regularization parameter in LASSO.
Moreover, when considering the secure implementation of cost-effective
computers, reducing the computation time is essential to achieve a resource-
aware encrypted controller by streamlining the proposed controller.

## References

  * [1] Ching-Ping Chou and B. Hannaford, “Measurement and modeling of Mckibben pneumatic artificial muscles,” _IEEE Transactions on Robotics and Automation_ , vol. 12, no. 1, pp. 90–102, 1996.
  * [2] D. G. Caldwell, G. A. Medrano-Cerda, and M. Goodwin, “Control of pneumatic muscle actuators,” _IEEE Control Systems Magazine_ , vol. 15, no. 1, pp. 40–48, 1995.
  * [3] B. Tondu and P. Lopez, “Modeling and control of mckibben artificial muscle robot actuators,” _IEEE Control Systems Magazine_ , vol. 20, no. 2, pp. 15–38, 2000.
  * [4] T. V. Minh, B. Kamers, T. Tjahjowidodo, H. Ramon, and H. Van Brussel, “Modeling torque-angle hysteresis in a pneumatic muscle manipulator,” in _IEEE/ASME International Conference on Advanced Intelligent Mechatronics_ , 2010, pp. 1122–1127.
  * [5] P. Beyl, K. Knaepen, S. Duerinck, M. V. Damme, B. Vanderborght, R. Meeusen, and D. Lefeber, “Safe and compliant guidance by a powered knee exoskeleton for robot-assisted rehabilitation of gait,” _Advanced Robotics_ , vol. 25, no. 5, pp. 513–535, 2011.
  * [6] G. Andrikopoulos, G. Nikolakopoulos, and S. Manesis, “Non-linear control of pneumatic artificial muscles,” in _Mediterranean Conference on Control and Automation_ , 2013, pp. 729–734.
  * [7] M. Tóthová and J. Pitel, “Simulation of actuator dynamics based on geometric model of pneumatic artificial muscle,” in _International Symposium on Intelligent Systems and Informatics_ , 2013, pp. 233–237.
  * [8] G. Andrikopoulos, G. Nikolakopoulos, and S. Manesis, “Advanced nonlinear PID-based antagonistic control for pneumatic muscle actuators,” _IEEE Transactions on Industrial Electronics_ , vol. 61, no. 12, pp. 6926–6937, 2014\.
  * [9] G. Andrikopoulos, G. Nikolakopoulos, and S. Manesis, “Novel considerations on static force modeling of pneumatic muscle actuators,” _IEEE/ASME Transactions on Mechatronics_ , vol. 21, no. 6, pp. 2647–2659, 2016.
  * [10] P. K. Jamwal, S. Hussain, Y. H. Tsoi, and S. Q. Xie, “Musculoskeletal model for path generation and modification of an ankle rehabilitation robot,” _IEEE Transactions on Human-Machine Systems_ , vol. 50, no. 5, pp. 373–383, 2020.
  * [11] D. Gong, L. Hao, J. Yu, and G. Zuo, “Bionic design of a dexterous anthropomorphic hand actuated by antagonistic PAMs,” in _IEEE International Conference on Real-time Computing and Robotics_ , 2020, pp. 493–498.
  * [12] K. Zhou, X. Sun, R. Yu, X. Liu, and B. Wang, “Pneumatic artificial muscle antagonistic joint trajectory tracking using adaptive explicit model predictive control,” in _International Conference on Computers and Automation_ , 2022, pp. 48–55.
  * [13] G. Andrikopoulos, G. Nikolakopoulos, and S. Manesis, “A survey on applications of pneumatic artificial muscles,” in _Mediterranean Conference on Control Automation_ , 2011, pp. 1439–1446.
  * [14] I. Sardellitti, G. Palli, N. G. Tsagarakis, and D. G. Caldwell, “Antagonistically actuated compliant joint: Torque and stiffness control,” in _IEEE/RSJ International Conference on Intelligent Robots and Systems_ , 2010, pp. 1909–1914.
  * [15] L. M. Sui and S. Q. Xie, “A model of pneumatic muscle actuated joint using linearized method,” in _International Conference on Mechatronics and Machine Vision in Practice_ , 2012, pp. 414–419.
  * [16] B. Ugurlu, P. Forni, C. Doppmann, and J. Morimoto, “Torque and variable stiffness control for antagonistically driven pneumatic muscle actuators via a stable force feedback controller,” in _IEEE/RSJ International Conference on Intelligent Robots and Systems_ , 2015, pp. 1633–1639.
  * [17] Y. Yamamoto, N. Matsunaga, and H. Okajima, “Robust variable stiffness control of McKibben type pneumatic artificial muscle arm by using multiple model error compensators,” in _International Conference on Control, Automation and Systems_ , 2017, pp. 957–962.
  * [18] X. Zhao, H. Ma, D. Ye, and D. Zhang, “Independent stiffness and force control of antagonistic pneumatic artificial muscles joint,” in _IEEE International Conference on Advanced Robotics and Mechatronics_ , 2017, pp. 734–739.
  * [19] D. Tanaka, H. Maeda, and T. Nakamura, “Joint stiffness and position control of an artificial muscle manipulator considering instantaneous load,” in _Annual Conference of IEEE Industrial Electronics_ , 2009, pp. 2259–2264.
  * [20] T. Choi, B. Choi, and K. Seo, “Position and compliance control of a pneumatic muscle actuated manipulator for enhanced safety,” _IEEE Transactions on Control Systems Technology_ , vol. 19, no. 4, pp. 832–842, 2011\.
  * [21] J. Cao, S. Q. Xie, and R. Das, “MIMO sliding mode controller for gait exoskeleton driven by pneumatic muscles,” _IEEE Transactions on Control Systems Technology_ , vol. 26, no. 1, pp. 274–281, 2018.
  * [22] N. Saito and T. Satoh, “Posture control considering joint stiffness of a robotic arm driven by rubberless artificial muscle,” _International Journal of Automation Technology_ , vol. 10, no. 4, pp. 503–510, 2016.
  * [23] B. Ugurlu, P. Forni, C. Doppmann, E. Sariyildiz, and J. Morimoto, “Stable control of force, position, and stiffness for robot joints powered via pneumatic muscles,” _IEEE Transactions on Industrial Informatics_ , vol. 15, no. 12, pp. 6270–6279, 2019.
  * [24] T. Shin and K. Kogiso, “Sensorless angle and stiffness control of antagonistic PAM actuator using reference set,” _Advanced Robotics_ , vol. 36, no. 9, pp. 423–437, 2022.
  * [25] H. Li, K. Tadano, and K. Kawashima, “Achieving force perception in master-slave manipulators using pneumatic artificial muscles,” in _SICE Annual Conference_ , 2012, pp. 1342–1345.
  * [26] T. Sasaki, T. Nagai, and K. Kawashima, “Remote control of backhoe for rescue activities using pneumatic robot system,” in _IEEE International Conference on Robotics and Automation_ , 2006, pp. 3177–3182.
  * [27] T. Kato, T. Higashi, and K. Shimizu, “Teleoperation of a robot arm system using pneumatic artificial rubber muscles: Teleoperation over the internet using udp and a web camera,” in _International Conference on Broadband, Wireless Computing, Communication and Applications_ , 2010, pp. 714–718.
  * [28] H. Toda, M. Tada, T. Maruyama, and Y. Kurita, “Effect of contraction parameters on swing support during walking using wireless pneumatic artificial muscle driver: A preliminary study,” in _Annual Conference of the Society of Instrument and Control Engineers_ , 2019, pp. 727–732.
  * [29] S. Yean, M. Tada, H. Toda, B.-S. Lee, and Y. Kurita, “Adaptive automatic controller for swing assist by pneumatic artificial muscle,” in _IEEE Sensors Applications Symposium_ , 2020, pp. 1–6.
  * [30] R. Langner, “Stuxnet: Dissecting a cyberwarfare weapon,” _IEEE Security & Privacy_, vol. 9, no. 3, pp. 49–51, 2011.
  * [31] S. Gorman, Y. J. Dreazen, and Cole, “Insurgents hack US drones,” https://www.wsj.com/articles/SB126102247889095011, 2009.
  * [32] D. Quarta, M. Pogliani, M. Polino, F. Maggi, A. M. Zanchettin, and S. Zanero, “An experimental security analysis of an industrial robot controller,” in _IEEE Symposium on Security and Privacy_ , 2017, pp. 268–286.
  * [33] K. Kogiso and T. Fujita, “Cyber-security enhancement of networked control systems using homomorphic encryption,” in _IEEE Conference on Decision and Control_ , 2015, pp. 6836–6843.
  * [34] K. Teranishi, K. Kogiso, and J. Ueda, “Encrypted feedback linearization and motion control for manipulator with somewhat homomorphic encryption,” in _IEEE/ASME International Conference on Advanced Intelligent Mechatronics_ , 2020, pp. 613–618.
  * [35] Y. Qiu and J. Ueda, “Encrypted motion control of a teleoperation system with security-enhanced controller by deception,” in _Dynamic Systems and Control Conference_ , vol. 59148, 2019, p. V001T07A006.
  * [36] T. Shin, K. Teranishi, and K. Kogiso, “Cyber-secure pneumatic actuator system equipped with encrypted controller and attack detectors,” _Advanced Robotics_ , vol. 36, no. 9, pp. 438–449, 2022.
  * [37] R. Baba, K. Kogiso, and M. Kishida, “Detection method of controller falsification attacks against encrypted control system,” in _SICE Annual Conference_ , 2018, pp. 244–248.
  * [38] M. Miyamoto, K. Teranishi, K. Emura, and K. Kogiso, “Cybersecurity-enhanced encrypted control system using keyed-homomorphic public key encryption,” _IEEE Access_ , vol. 11, pp. 45 749–45 760, 2023.
  * [39] M. Schulze Darup, “Encrypted polynomial control based on tailored two-party computation,” _International Journal of Robust and Nonlinear Control_ , vol. 30, no. 11, pp. 4168–4187, 2020.
  * [40] T. Shin, T. Ibayashi, and K. Kogiso, “Detailed dynamic model of antagonistic PAM system and its experimental validation: Sensorless angle and torque control with UKF,” _IEEE/ASME Transactions on Mechatronics_ , vol. 27, no. 3, pp. 1715–1726, 2022.
  * [41] Y. Takeda and K. Kogiso, “Experimental regularization parameter search for polynomial approximation of nonlinear PAM controller,” in _IEEE/SICE International Symposium on System Integration_ , 2023, pp. 963–968.

|  Yuta Takeda received the B.S. degree in Informatics and Engineering from
The University of Electro-Communications, Tokyo, Japan, in 2022. He is
currently an M.S. student at The University of Electro-Communications, Tokyo,
Japan. His research interests include control applications and
modeling/control of pneumatic artificial muscles.
---|---
|  Takaya Shin received the B.S. and M.S. degrees in Informatics and
Engineering from The University of Electro-Communications, Tokyo, Japan, in
2020 and 2022, respectively. He joined Daihen Co., Osaka, Japan. His research
interests include control applications and modeling/control of pneumatic
artificial muscles.
---|---
|  Kaoru Teranishi received the B.S. degree in electromechanical engineering
from National Institute of Technology, Ishikawa College, Ishikawa, Japan, in
2019. He also obtained the M.S. degree in Mechanical and Intelligent Systems
Engineering from The University of Electro-Communications, Tokyo, Japan, in
2021. He is currently a Ph.D. student at The University of Electro-
Communications. From October 2019 to September 2020, he was a visiting scholar
of the Georgia Institute of Technology, GA, USA. Since April 2021, he has been
a Research Fellow of the Japan Society for the Promotion of Science. His
research interests include control theory and cryptography for cyber-security
of control systems.
---|---
|  Kiminao Kogiso received the B.E., M.E., and Ph.D. degrees in mechanical
engineering from Osaka University, Japan, in 1999, 2001, and 2004,
respectively. He was appointed a postdoctoral fellow in the 21st Century COE
Program and as an Assistant Professor in the Graduate School of Information
Science, Nara Institute of Science and Technology, Nara, Japan, in April 2004
and July 2005, respectively. From November 2010 to December 2011, he was a
visiting scholar at the Georgia Institute of Technology, GA, USA. In March
2014, he was promoted to associate professor in the Department of Mechanical
and Intelligent Systems Engineering at The University of Electro-
Communications, Tokyo, Japan. Since April 2023, he has been a professor in the
same department. His research interests include the cybersecurity of control
systems, constrained control, control of decision-makers, and their
applications.
---|---


3.5cm3.5cm1.3cm2cm

# On the containment problem and sporadic simplicial line arrangements

Marek Janasz

###### Abstract

In the paper we present two examples of inductively free sporadic simplicial
arrangements of $31$ lines that are non-isomorphic, which allow us to answer
negatively questions on the containment problem recently formulated by Drabkin
and Seceleanu.

Keywords simplicial arrangements, homogeneous ideals, symbolic powers

Mathematics Subject Classification (2000) 14N20, 13A15, 52A20

##  1 Introduction

Study of the relation between symbolic and ordinary power of homogeneous
ideals in the polynomial ring over a given field $\mathbb{K}$ has a long
history and is derived from many different problems in mathematics. In $1995$,
Eisenbud and Mazur in [17], referring to the proof of Fermat’s Last Theorem,
were investigating the so-called ”fitting ideals” and some symbolic powers of
certain associated ideals. What they proved, among other things, is that
$I^{(2)}\subset\mathfrak{m}I$ in the case of perfect ideals of codimension
$2$. They also showed that this kind of the containment holds for several
other classes of ideals. In $2013$, Harbourne and Huneke in [23] proposed a
certain generalization and began to study the relation
$I^{(m)}\subset\mathfrak{m}^{k}I^{r}$, and their work was continued in [2, 7].

Another way, which in fact ends in the investigations on symbolic and ordinary
powers of ideals, were the articles by Skoda [29] and Waldschmidt [33]. They
focused on some estimates of the degree of hypersurfaces in
$\mathbb{P}^{N}_{\mathbb{K}}$ passing through fixed points with prescribed
multiplicities. A paper by Chudnovski [5] fits into these considerations.
Using some complex analysis tools, Chudnovski improved the results of Skoda
and Waldschmidt in $\mathbb{P}^{2}$ and he formulated a still open conjecture
for the case of $\mathbb{P}^{N}$. The generalization of this conjecture was
also given by Demailly in [10], and the combination of these conjectures is
the subject of intense research [2, 3, 14, 16, 20, 27].

Using the Nagata-Zariski theorem makes it possible to relate geometric
questions to algebra. Therefore, the study of containment relations between
ordinary and symbolic powers of homogeneous ideals of points is a connection
between conjectures formulated by Chudnovsky and Demailly. This perspective
led to an increased interest in the so-called containment problem, i.e., the
determination of the exponents $(m,r)$ for which the $m$-th symbolic power of
a homogeneous ideal is contained in the $r$-th ordinary power of that ideal.
An initial work on this topic was begun by Hochster in $1973$ [24], but the
groundbreaking result was published only in $2001$ by Ein, Lazarsfeld, and
Smith [18], where they gave a lower bound on the exponent of the symbolic
power, which depends on the dimension of the space $\mathbb{P}^{N}$. Since
then, the cases unsolved by the aforementioned theorem have become the subject
of intensive study, in particular the smallest case from the perspective of
the magnitude of the powers, namely the containment $I^{(3)}\subset I^{2}$ for
the ideals of reduced points in $\mathbb{P}^{2}$. While at the beginning the
researchers tried to prove that this particular containment holds for all
homogeneous ideals, after the paper [15], where the first counterexample
defined over $\mathbb{C}$ has been presented, a lot of counterexamples defined
over different fields have been published (see [4, 26]). Despite a growing
number of counterexamples, the true nature of the relation between $I^{(3)}$
and $I^{2}$ is still unknown. In [11], Drabkin and Seceleanu study reflection
arrangements given by (irreducible) complex pseudoreflection groups. As a
result, they give a complete description of the relation between the third
symbolic power and the second ordinary power of radical ideal
$J({\mathcal{A}})$, which defines the singular locus of the complex reflection
arrangement ${\mathcal{A}}$. The work on this problem motivates them to state
some open questions, among which they ask ([11, Question 6.7.-6.8.]): Are the
containments $(J({\mathcal{A}}))^{(2r-1)}\subseteq(J({\mathcal{A}}))^{r}$
always satisfied for any $r\geq 2$ and any hyperplane arrangement that is
inductively/recursively free?

In the present paper we give a negative answer to these questions, namely we
prove the following.

######  Theorem A.

There are two non-isomorphic inductively free simplicial arrangements
consisting of $31$ lines such that they have the same weak combinatorics, and
having the property that for one arrangement the containment
$(J({\mathcal{A}}))^{(3)}\subseteq(J({\mathcal{A}}))^{2}$ holds, but does not
for the other.

The structure of the paper is as follows. In Section 2, we recall some basic
definitions and tools that we will use in the rest of this paper concerning
line arrangements and symbolic powers of homogeneous ideals. In Section 3 we
give very detailed information about a family of line arrangement known as
${\mathcal{A}}(12k+7)$, giving line equations and proving that some line
arrangements from this family are inductively free. This result is used in
Section 4, where we prove Main Theorem of this paper. At the end, we provide
our `SINGULAR` code to let interested readers check the containment between
$(J({\mathcal{A}}))^{(3)}$ and $(J({\mathcal{A}}))^{2}$.

##  2 Preliminaries

In this section we recall all necessary definitions regarding hyperplane
arrangements that we will exploit in the paper. For more information regarding
this subject, please consult [13, 28].

Let $\mathbb{K}$ be a field of characteristic zero and let $V$ be a fixed
vector space of dimension $\ell$ over $\mathbb{K}$. Let
$\\{x_{1},\ldots,x_{\ell}\\}$ be the dual basis of $V^{*}$ associated with
$V$, then the symmetric algebra $S(V^{*})$ is isomorphic to the ring of
polynomials $S=\mathbb{K}[x_{1},\ldots,x_{\ell}].$

A pair $({\mathcal{A}},V)$ is called an $\ell$-arrangement of hyperplanes,
i.e., this is an arrangement of dimension $(\ell-1)$ linear subspaces in $V$.
The symbol $\Phi_{\ell}$ denotes the empty $\ell$-arrangement. If the
dimension is clear from the context, we use the name arrangement for short.
Each hyperplane $H\in{\mathcal{A}}$ is the kernel (up to a constant) of a
linear form $l_{H}\in V^{*}$. The product of all linear forms

$Q({\mathcal{A}})=\prod_{H\in{\mathcal{A}}}l_{H}$

is called the defining polynomial of ${\mathcal{A}}$. In the case of empty
arrangements, we put $Q(\Phi_{l})=1$.

By $L({\mathcal{A}})$ we denote the intersection lattice of ${\mathcal{A}}$,
i.e., the set of all non-empty intersections of hyperplanes $H_{i}$ in
${\mathcal{A}}$. Taking any $X\in L({\mathcal{A}})$, a subarrangement
${\mathcal{A}}_{X}$ of ${\mathcal{A}}$ is called localization and it is
defined as

${\mathcal{A}}_{X}=\\{H\in{\mathcal{A}}\>|\>X\subseteq H\\}.$

For a chosen $X$, we define a subarrangement of ${\mathcal{A}}$ by

${\mathcal{A}}^{X}=\\{X\cap H\>:X\not\subseteq H\>\text{and}\>X\cap
H\not=\emptyset\\},$

which we call the restriction of ${\mathcal{A}}$ to $X$.

######  Definition 2.1.

A simplicial arrangement is a finite set
${\mathcal{A}}=\\{H_{1},\ldots,H_{n}\\}$ of (central) hyperplanes in
$\mathbb{R}^{\ell}$ such that all connected components of the complement

$M(\mathcal{A}):=\mathbb{R}^{\ell}\setminus\bigcup_{H\in{\mathcal{A}}}H$

are simplicial cones.

Denote by $\mathop{Der}\nolimits_{\mathbb{K}}(S)$ the set of all
$\mathbb{K}$-linear maps (derivations) $\theta:S\longrightarrow S$ such that
for all $f,g\in S$ one has

$\theta(fg)=f\theta(g)+g\theta(f).$

It is known that the set $\left\\{\tfrac{\partial}{\partial
x_{i}}\right\\}_{i=1}^{\ell}$ forms a (canonical) basis for
$\mathop{Der}\nolimits_{\mathbb{K}}(S)$, i.e.,

$\mathop{Der}\nolimits_{\mathbb{K}}(S)=\bigoplus_{i=1}^{\ell}S\cdot\tfrac{\partial}{\partial
x_{i}}.$

Any homogeneous element $0\not=\theta\in\mathop{Der}\nolimits_{\mathbb{K}}(S)$
can be expressed as
$\theta=\sum_{i=1}^{\ell}g_{i}\cdot\tfrac{\partial}{\partial x_{i}}$, where
$g_{i}\in S$ are homogeneous polynomials of degree $d$. For such $\theta$ we
denote by $\mathop{pdeg}\nolimits\theta=d$ its polynomial degree.

For any $f\in S$ being homogeneous, we define an $S$-submodule of
$\mathop{Der}\nolimits_{\mathbb{K}}(S)$ as

$D(f)=\\{\theta\in\mathop{Der}\nolimits_{\mathbb{K}}(S):\theta(f)\in f\cdot
S\\}.$

In the case of arrangement ${\mathcal{A}}$, we use the notation
$D({\mathcal{A}})$ instead of $D(Q({\mathcal{A}}))$.

######  Definition 2.2.

If $D({\mathcal{A}})$ is a free $S$-module, then we say that ${\mathcal{A}}$
is a _free arrangement_.

Let ${\mathcal{A}}$ be a free arrangement for which
$\\{\theta_{1},\ldots,\theta_{\ell}\\}$ is a homogeneous basis of
$D({\mathcal{A}})$. We say that the set of integers
$\exp({\mathcal{A}})=\\{\mathop{pdeg}\nolimits\theta_{1},\ldots,\mathop{pdeg}\nolimits\theta_{\ell}\\}$
is the set of the exponents of ${\mathcal{A}}$.

If we denote by $\theta_{E}\in\mathop{Der}\nolimits_{\mathbb{K}}(S)$ the Euler
derivation, then we have the decomposition of $D({\mathcal{A}})$, namely

$D({\mathcal{A}})=S\cdot\theta_{E}\oplus D_{0}({\mathcal{A}}).$

Now for an arrangement ${\mathcal{A}}$ and fixed $H\in{\mathcal{A}}$, it is
convenient to study triples of arrangements
$({\mathcal{A}},{\mathcal{A}}^{{}^{\prime}},{\mathcal{A}}^{{}^{\prime\prime}})$
of arrangements, where
${\mathcal{A}}^{{}^{\prime}}={\mathcal{A}}\setminus\\{H\\}$ and
${\mathcal{A}}^{{}^{\prime\prime}}={\mathcal{A}}^{H}.$ The next theorem is
very useful in all our considerations.

######  Theorem 2.3.

(Addition-Deletion, see [28]) Suppose $\mathcal{A}\not=\Phi_{\ell}$ . Let
$({\mathcal{A}},{\mathcal{A}}^{{}^{\prime}},{\mathcal{A}}^{{}^{\prime\prime}})$
be a triple. Any two of the following statements imply the third:

$\displaystyle\mathcal{A}\text{ is free with }{\rm exp}(\mathcal{A})$
$\displaystyle=\\{b_{1},\ldots,b_{\ell-1},b_{\ell}\\},$
$\displaystyle\mathcal{A}^{{}^{\prime}}\text{ is free with }{\rm
exp}(\mathcal{A}^{{}^{\prime}})$
$\displaystyle=\\{b_{1},\ldots,b_{\ell-1},b_{\ell}-1\\},$
$\displaystyle\mathcal{A}^{{}^{\prime\prime}}\text{ is free with }{\rm
exp}(\mathcal{A}^{{}^{\prime\prime}})$
$\displaystyle=\\{b_{1},\ldots.,b_{\ell-1}\\}.$

In this paper we deal with line arrangements ${\mathcal{A}}$ defined over the
complex numbers, therefore we will use the following reformulation of Theorem
2.3.

######  Theorem 2.4.

Let ${\mathcal{A}}$ be a line arrangement in $\mathbb{P}^{2}_{\mathbb{C}}$ and
$H\in{\mathcal{A}}$. Let
${\mathcal{A}}^{\prime}:={\mathcal{A}}\setminus\\{H\\}$. If the following
conditions hold:

  1. (1)

${\mathcal{A}}^{\prime}$ is free and has the exponents ${\rm
exp}(\mathcal{A}^{\prime})=\\{1,a,b\\}$,

  2. (2)

$|Sing({\mathcal{A}})\cap H|=b+1$ (or $a+1$, respectively),

then ${\mathcal{A}}$ is free with the exponents ${\rm
exp}({\mathcal{A}})=\\{1,a+1,b\\}$ (or ${\rm
exp}({\mathcal{A}})=\\{1,a,b+1\\}$, respectively).

In the sequel, we focus on the following definition.

######  Definition 2.5.

([28, Definition 4.53]). The class $\mathcal{IF}$ of inductively free
arrangements is the smallest class of arrangements which satisfies both
conditions:

  1. (1)

$\Phi_{\ell}\in\mathcal{IF}$ for $\ell\geq 0$,

  2. (2)

if there exists $H\in{\mathcal{A}}$ such that
${\mathcal{A}}^{\prime\prime}\in\mathcal{IF}$,
${\mathcal{A}}^{\prime}\in\mathcal{IF}$, and ${\rm
exp}({\mathcal{A}}^{\prime\prime})\subset{\rm exp}({\mathcal{A}}^{\prime})$,
then ${\mathcal{A}}\in\mathcal{IF}$.

##  3 Examples of inductively free arrangements

The main object of our considerations is a special family of line
arrangements, denoted in the literature by ${\mathcal{A}}(12k+7)$. This
infinite family was originally described in the paper by Grünbaum [1]. Here,
we recall this construction and its basic properties.

For fixed $k$, each element of the family consists of exactly $12k+7$ lines,
including the line at infinity $z=0$. The equations of these lines are given
explicitly in Table 1.

${\mathcal{A}}(12k+7)$
---
$2x-eiz$, |
$x-ey+iez$, | for $i\in\\{-(k+1),-k,\ldots,-1,0,1,\ldots,k,k+1\\}$
$x+ey-iez$, |
$2y-jz$, |
$ex-y+jz$, | for $j\in\\{-(k-1),-(k-2),\ldots,-1,0,1,\ldots,k-2,k-1\\}$
$ex+y-jz$, |
$z$ |
Table 1:  Equations of lines of ${\mathcal{A}}(12k+7)$.

The arrangements ${\mathcal{A}}(19)$ and ${\mathcal{A}}(31)$ are exactly the
sporadic simplicial arrangements $A(19,1)$ and $A(31,2)$ listed in Grünbaum’s
catalogue [22]. From [12], we know that these are free arrangements with the
exponents ${\rm exp}({\mathcal{A}}(19))=\\{1,7,11\\}$ and ${\rm
exp}({\mathcal{A}}(31))=\\{1,13,17\\}$. We start with the following
combinatorial observation that is crucial for our further considerations.
According to our best knowledge, this observation is not known in the
literature.

######  Theorem 3.1.

The arrangements of line ${\mathcal{A}}(19,1)$ and ${\mathcal{A}}(31,2)$ are
inductively free.

###### Proof.

We will divide our proof of this theorem into two steps. In the first step, we
will show that the arrangement ${\mathcal{A}}(19,1)$ is inductively free. For
this purpose, we present Table 2 below, where we deliver the sequences of the
exponents for ${\mathcal{A}}^{{}^{\prime}}$, the equation of each line that we
add to the arrangement starting from $\Phi_{3}$, and then the exponents of
${\rm exp}({\mathcal{A}}^{{}^{\prime\prime}})$. Each subsequent row of the
table allows us to verify the conditions contained in Theorem 2.4.

$\exp\>{\mathcal{A}}^{{}^{\prime}}$ | $\ell_{i}$ | $\exp\>{\mathcal{A}}^{{}^{\prime\prime}}$ | $\exp\>{\mathcal{A}}^{{}^{\prime}}$ | $\ell_{i}$ | $\exp\>{\mathcal{A}}^{{}^{\prime\prime}}$
---|---|---|---|---|---
$\\{0,0,0\\}$ | $\Phi_{3}$ | $\\{0,0\\}$ | $\\{1,4,5\\}$ | $\ell_{10}:x+ey+ez$ | $\\{1,5\\}$
$\\{0,0,1\\}$ | $\ell_{1}:z$ | $\\{0,1\\}$ | $\\{1,5,5\\}$ | $\ell_{11}:x-ey-ez$ | $\\{1,5\\}$
$\\{0,1,1\\}$ | $\ell_{2}:ex+y$ | $\\{1,1\\}$ | $\\{1,5,6\\}$ | $\ell_{12}:x+ey-ez$ | $\\{1,5\\}$
$\\{1,1,1\\}$ | $\ell_{3}:ex-y$ | $\\{1,1\\}$ | $\\{1,5,7\\}$ | $\ell_{13}:x-ey+ez$ | $\\{1,7\\}$
$\\{1,1,2\\}$ | $\ell_{4}:y$ | $\\{1,1\\}$ | $\\{1,6,7\\}$ | $\ell_{14}:2x-2ez$ | $\\{1,7\\}$
$\\{1,1,3\\}$ | $\ell_{5}:x$ | $\\{1,1\\}$ | $\\{1,7,7\\}$ | $\ell_{15}:x+ez$ | $\\{1,7\\}$
$\\{1,1,4\\}$ | $\ell_{6}:x+ey$ | $\\{1,1\\}$ | $\\{1,7,8\\}$ | $\ell_{16}:x+ey+2ez$ | $\\{1,7\\}$
$\\{1,1,5\\}$ | $\ell_{7}:x-ey$ | $\\{1,5\\}$ | $\\{1,7,9\\}$ | $\ell_{17}:x-ey+2ez$ | $\\{1,7\\}$
$\\{1,2,5\\}$ | $\ell_{8}:2x-ez$ | $\\{1,5\\}$ | $\\{1,7,10\\}$ | $\ell_{18}:x-ey-2ez$ | $\\{1,7\\}$
$\\{1,3,5\\}$ | $\ell_{9}:2x+ez$ | $\\{1,5\\}$ | $\\{1,7,11\\}$ | $\ell_{19}:x+ey-2ez$ | $\\{1,7\\}$
Table 2: . List of the exponents for arrangements building
${\mathcal{A}}(19)$, where $e=\sqrt{3}$.

Based on the last row in Table 2, we see that the arrangement
${\mathcal{A}}(19)$ is a free with the exponents ${\rm
exp}(\mathcal{A}(19))=\\{1,7,11\\}$.

For the second part of the proof, we apply Theorem 2.4 to arrangement
$\mathcal{A}(19)$ by adding suitably chosen lines that are indicated in Table
3 below.

$\ell_{19+i}$ | Equations of lines | $\ell_{19+i}$ | Equations of lines
---|---|---|---
$\ell_{20}:$ | $ex+y-z$ | $\ell_{26}:$ | $x+ey+3ez$
$\ell_{21}:$ | $ex+y+z$ | $\ell_{27}:$ | $x+ey-3ez$
$\ell_{22}:$ | $ex-y+z$ | $\ell_{28}:$ | $x-ey+3ez$
$\ell_{23}:$ | $ex-y-z$ | $\ell_{29}:$ | $x-ey-3ez$
$\ell_{24}:$ | $2y-z$ | $\ell_{30}:$ | $2x+3ez$
$\ell_{25}:$ | $2y+z$ | $\ell_{31}:$ | $2x-3ez$
Table 3: . Equations of lines $\ell_{19+i}$ with $i\in\\{1,...,12\\}$.

Attaching successively lines $\ell_{20},...,\ell_{31}$ to $\mathcal{A}(19)$,
and using each time Theorem 2.4, we conclude that the obtained arrangements
are free. All the details describing our procedure are presented in the
following diagram below:

$\begin{array}[]{c}{\mathcal{A}}(19)\\\ {\rm
exp}({\mathcal{A}})=\\{1,6,11\\}\end{array}\longrightarrow\begin{array}[]{c}{\mathcal{A}}(19)\cup\\{\ell_{20}\\}\\\
{\rm exp}({\mathcal{A}})=\\{1,8,11\\}\end{array}\longrightarrow\ldots$
$\longrightarrow\begin{array}[]{c}{\mathcal{A}}(19)\cup\\{\ell_{20},\ldots,\ell_{24}\\}\\\
{\rm
exp}({\mathcal{A}})=\\{1,11,11\\}\end{array}\longrightarrow\begin{array}[]{c}{\mathcal{A}}(19)\cup\\{\ell_{20},\ldots,\ell_{25}\\}\\\
{\rm exp}({\mathcal{A}})=\\{1,11,12\\}\end{array}\longrightarrow$
$\longrightarrow\begin{array}[]{c}{\mathcal{A}}(19)\cup\\{\ell_{20},\ldots,\ell_{26}\\}\\\
{\rm
exp}({\mathcal{A}})=\\{1,11,13\\}\end{array}\longrightarrow\begin{array}[]{c}{\mathcal{A}}(19)\cup\\{\ell_{20},\ldots,\ell_{27}\\}\\\
{\rm exp}({\mathcal{A}})=\\{1,12,13\\}\end{array}$

$\longrightarrow\begin{array}[]{c}{\mathcal{A}}(19)\cup\\{\ell_{20},\ldots,\ell_{28}\\}\\\
{\rm
exp}({\mathcal{A}})=\\{1,13,13\\}\end{array}\longrightarrow\begin{array}[]{c}{\mathcal{A}}(19)\cup\\{\ell_{20},\ldots,\ell_{29}\\}\\\
{\rm exp}({\mathcal{A}})=\\{1,13,14\\}\end{array}\longrightarrow$
$\ldots\longrightarrow\begin{array}[]{c}{\mathcal{A}}(19)\cup\\{\ell_{20},\ldots,\ell_{31}\\}\\\
{\rm exp}({\mathcal{A}})=\\{1,13,17\\}\end{array}.$

Thus we obtain that the arrangement

${\mathcal{A}}(31):={\mathcal{A}}(19)\cup\\{\ell_{20},\ldots,\ell_{31}\\}$

is inductively free with the exponents ${\rm
exp}({\mathcal{A}}(31))=\\{1,13,17\\},$ which completes the proof. ∎

##  4 Inductively free arrangements and a counterexample to the containment
problem.

The line arrangement ${\mathcal{A}}(31)$ from the previous part of our paper
turns out to be very important in the context of an open problem in the so-
called containment problems for symbolic powers of homogeneous ideals. Let us
recall here the definition of symbolic powers.

######  Definition 4.1.

Let $I\subseteq S$ be a homogeneous ideal. For a fixed positive integer $m$,
we define the $m$-th symbolic power of the ideal $I$ as

$I^{(m)}=S\cap\Big{(}\bigcap_{Q\in\text{Ass}(I)}I^{m}_{Q}\Big{)},$

where $Ass(I)$ denotes the set of all prime ideals associated with I and
$I_{Q}$ denotes the location of $I$ at $Q$.

Reader unfamiliar with symbolic power and containment problem is referred to
[30].

In [11], Drabkin and Seceleanu study arrangements of hyperplanes that come
from irreducible complex reflection groups, proving that for some cases we
have the failure of the containment

$(J({\mathcal{A}}))^{(3)}\not\subseteq(J({\mathcal{A}}))^{2},$ (1)

where $J(\mathcal{A})$ denotes the radical ideal associated with the
configuration of all intersection points of a given arrangement $\mathcal{A}$
and $(J({\mathcal{A}}))^{(3)}$ denotes the third symbolic power of
$J({\mathcal{A}})$. In the light of the results obtained in [11], it is
natural to ask the following question.

######  Question 4.2.

([11, Question 6.7.]) Are the containments
$(J({\mathcal{A}}))^{(2r-1)}\subseteq(J({\mathcal{A}}))^{r}$ always satisfied
for any $r\geq 2$ and any hyperplane arrangement that is inductively free?

Here we answer negative to this question for $r=2$ and by taking the whole
singular locus of the arrangement $\mathcal{A}(31)$ described in the previous
section. In fact, we will show even more, namely we will present two
inductively free arrangements of $31$ lines such that for one arrangement the
condition (1) holds, but for the second one it does not. In order to do so,
let us present briefly a construction of the second arrangement of $31$ lines.
Surprisingly, this is a simplicial line arrangement which is denoted by
${\mathcal{A}}(31,3)$ in Grünbaum’s catalogue [22]. Here are the details.

Consider ten lines given by:

$x\pm aez,\>2y\pm az,\>2x\pm bz,$

where by $e=\sqrt{3}$, $a\in\\{0,1\\}$, and $b\in\\{0,1,3\\}$. The
visualization in the affine part of the projective plane of these ten lines is
presented in Figure 1.

$x=0$$2x=-ez$$2x=ez$$x=-ez$$x=ez$$2x=-3ez$$2x=3ez$$2y=0$$2y=z$$2y=-z$ Figure
1: . A set of $10$ lines inducing realizations of ${\mathcal{A}}(31,3)$.

Now we perform two rotation of these lines (affinely), firstly by angle
$60^{\circ}$, then by $120^{\circ}$, around the point $(0,0)$. In this way we
obtain $30$ lines. Finally, to obtain $31$ lines, we take to our arrangement
the line at infinity $z=0$. The resulting configuration is shown in Figure 2.

$\infty$ Figure 2: . An affine realization of the arrangement
${\mathcal{A}}(31,3)$. The symbol $\infty$ denotes the line at infinity $z=0$.

The following statement shows that the answer to Question 4.2 is negative.

######  Theorem 4.3.

The arrangement ${\mathcal{A}}(31,3)$ is inductively free and one has
$(J({\mathcal{A}}(31,3)))^{(3)}\not\subseteq(J({\mathcal{A}}(31,3)))^{2}$.

###### Proof.

To prove that the configuration ${\mathcal{A}}(31,3)$ is inductively free, we
will create a table containing the exponents of ${\rm
exp}({\mathcal{A}}^{{}^{\prime}})$, the equations of lines that we add to the
arrangement ${\mathcal{A}}(19,1)$, and then the exponents of ${\rm
exp}({\mathcal{A}}^{{}^{\prime\prime}})$.

$\exp\>{\mathcal{A}}^{{}^{\prime}}$ | $\ell_{i}$ | $\exp\>{\mathcal{A}}^{{}^{\prime\prime}}$ | $\exp\>{\mathcal{A}}^{{}^{\prime}}$ | $\ell_{i}$ | $\exp\>{\mathcal{A}}^{{}^{\prime\prime}}$
---|---|---|---|---|---
$\\{1,8,11\\}$ | $\ell_{20}:ex+y+z$ | $\\{1,11\\}$ | $\\{1,12,13\\}$ | $\ell_{26}:4x+ez$ | $\\{1,13\\}$
$\\{1,9,11\\}$ | $\ell_{21}:ex+y-z$ | $\\{1,11\\}$ | $\\{1,13,13\\}$ | $\ell_{27}:4x-ez$ | $\\{1,13\\}$
$\\{1,10,11\\}$ | $\ell_{22}:ex-y+z$ | $\\{1,11\\}$ | $\\{1,13,14\\}$ | $\ell_{28}:2x-2ey+ez$ | $\\{1,13\\}$
$\\{1,11,11\\}$ | $\ell_{23}:ex-y-z$ | $\\{1,11\\}$ | $\\{1,13,15\\}$ | $\ell_{29}:2x-2ey-ez$ | $\\{1,13\\}$
$\\{1,11,12\\}$ | $\ell_{24}:2y-z$ | $\\{1,11\\}$ | $\\{1,13,16\\}$ | $\ell_{30}:2x+2ey-ez$ | $\\{1,13\\}$
$\\{1,11,13\\}$ | $\ell_{25}:2y+z$ | $\\{1,13\\}$ | $\\{1,13,17\\}$ | $\ell_{31}:2x+2ey+ez$ | $\\{1,13\\}$
Table 4: . Transition from ${\mathcal{A}}(19,1)$ to ${\mathcal{A}}(31,3)$.

Observe that each row in Table 4 allows us to verify condition (2) in
Definition 2.5. Evidence of noncontainment
$(J({\mathcal{A}}(31,3)))^{(3)}\not\subseteq(J({\mathcal{A}}(31,3)))^{2}$ has
been provided in [25] and due to this reason we refer to this paper for
details. The verification was performed using Singular. The affine part
($z=1$) of the element
$F\in(J({\mathcal{A}}(31,3)))^{(3)}\setminus(J({\mathcal{A}}(31,3)))^{2}$ is
shown in Figure 3. ∎

Figure 3: . The element
$F\in(J({\mathcal{A}}(31,3)))^{(3)}\setminus(J({\mathcal{A}}(31,3)))^{2}$
consists of $21$ lines and a curve of degree 12.

The next example illustrates the fact that being inductively free for an
arrangement ${\mathcal{A}}$ does not directly transfer into lack of
containment $(J({\mathcal{A}}))^{(3)}\subseteq(J({\mathcal{A}}))^{2}$.

######  Theorem 4.4.

There are two non-isomorphic inductively free simplicial arrangements
consisting of $31$ lines such that they have the same weak combinatorics, and
having the property that for one arrangement the containment
$(J({\mathcal{A}}))^{(3)}\subseteq(J({\mathcal{A}}))^{2}$ holds, not for the
other.

In other words, the weak combinatorics of line arrangements does not determine
the property of being an example for the non-containment. For the clarity of
the exposition, let us recall that for an arrangement $\mathcal{A}$ of $d$
lines in the plane, the weak combinatorics is the vector of the form
$(d;t_{2},...,t_{d})$ with $t_{i}$ being the number of $i$-fold intersection
points in $\mathcal{A}$.

It is worth noticing that Theorem 4.4 should be compared with a result from
[19], where the authors observed a similar phenomenon, but in a different
setting, namely in the case of real line arrangements possessing the maximal
possible number of triple intersection points.

Let us present our proof of Theorem 4.4.

###### Proof.

The arrangement ${\mathcal{A}}(31,3)$ is constructed from configuration
${\mathcal{A}}(19,1)$ by adding $12$ appropriately chosen lines, according to
Table 4. The arrangements ${\mathcal{A}}(31,2)$ and ${\mathcal{A}}(31,3)$ are
not isomorphic simplicial arrangements of lines, and this fact is proved in
[22], even though they have the same weak combinatorics, namely:

$t_{2}=54,\>t_{3}=42,\>t_{4}=21,\>t_{5}=6,\>t_{6}=1,\>t_{8}=3,$

and $t_{i}=0$ for others.

Let us pass to the containment question
$(J({\mathcal{A}}))^{(3)}\subseteq(J({\mathcal{A}}))^{2}$. In Theorem 4.3, we
explained the non-containment for the singular locus of $\mathcal{A}(31,3)$,
and this check was done using `SINGULAR`. In the case of $\mathcal{A}(31,2)$,
we can preform exactly the same computations showing that the containment

$(J({\mathcal{A}}(31,2))^{(3)}\subseteq(J({\mathcal{A}}(31,2))^{2}$

does hold, which completes the proof. ∎

In Appendix, you can find a script that can be run in `SINGULAR` which allows
us to verify the containment
$(J({\mathcal{A}}(31,2))^{(3)}\subseteq(J({\mathcal{A}}(31,2))^{2}$.

######  Remark 4.5.

Let us point our here that there is a way to extend the class of inductively
free arrangements $\mathcal{IF}$, namely we can add the following condition
([28, Definition 4.60 (3)]):

> if there exists $H\in{\mathcal{A}}$ such that
> ${\mathcal{A}}^{\prime\prime}\in\mathcal{IF}$,
> ${\mathcal{A}}\in\mathcal{IF}$, and ${\rm
> exp}({\mathcal{A}}^{\prime\prime})\subset{\rm exp}({\mathcal{A}})$, then
> ${\mathcal{A}}^{\prime}\in\mathcal{IF}$,

then we come to the class of recursively free hyperplanes arrangements
$\mathcal{RF}$.

It is known that we have the following relations (see [9, 31, 34])

$\text{inductively free }\subsetneq\text{ recursively free }\subsetneq\text{
free}.$

It turns out that our example of a pair of arrangements ${\mathcal{A}}(31,2)$
and ${\mathcal{A}}(31,3)$ allows us to answer a question posed by Drabkin and
Seceleanu in the negative ([11, Question 6.8]). More precisely, our example
shows that the containment $J(\mathcal{A})^{(3)}\subset J(\mathcal{A})^{2}$
does not hold for recursively free line arrangements $\mathcal{A}$. It is
still an open question if we can find a configuration of lines which is
recursively free, but not inductively free and gives negative answer.

### Acknowledgments.

I would like to thank Grzegorz Malara and Piotr Pokora for all their help,
valuable comments and inspiring discussions.

### Appendix.

    proc PtsIdeal(poly p, poly q, poly r) {
    Ψmatrix M[2][3]=p,q,r,x,y,z;
    Ψideal @I=minor(M,2);
    Ψreturn(std(@I));
    }
    option(redSB);
    ring R=(0,e),(x,y,z),dp;
    minpoly=e2-3;

    /* The list L contains the coordinates of the singular points
    of the arrangement A(31,2). */

    "loading arrangement A(31,2)...";
    list L=
    (-7/2e),-1/2,1,(7/2e),-1/2,1,(7/2e),1/2,1,(-7/2e),1/2,1,
    (-3/2e),-11/2,1,(2e),5,1,(3/2e),11/2,1,(-2e),-5,1,
    (3/2e),-11/2,1,(-2e),5,1,(-3/2e),11/2,1,(2e),-5,1,
    (3/4e),5/4,1,(1/4e),7/4,1,(1/4e),-7/4,1,(3/4e),-5/4,1,
    (-3/4e),5/4,1,(-1/4e),7/4,1,(-1/4e),-7/4,1,(-3/4e),-5/4,1,
    (-e),-1/2,1,(-e),1/2,1,(e),-1/2,1,(e),1/2,1,
    (5/2e),-1/2,1,(-3/2e),7/2,1,(5/2e),1/2,1,(-3/2e),-7/2,1,
    (-5/2e),-1/2,1,(3/2e),7/2,1,(-5/2e),1/2,1,(3/2e),-7/2,1,
    (-e),-4,1,(-e),4,1,(e),4,1,(e),-4,1,
    (3/2e),5/2,1,(-2e),-1,1,(-3/2e),5/2,1,(2e),-1,1,
    (-3/2e),-5/2,1,(2e),1,1,(3/2e),-5/2,1,(-2e),1,1,
    (-1/2e),-7/2,1,(-1/2e),7/2,1,(1/2e),7/2,1,(1/2e),-7/2,1,
    (-3/2e),-1/2,1,(e),2,1,(3/2e),-1/2,1,(-e),2,1,
    (3/2e),1/2,1,(-e),-2,1,(-3/2e),1/2,1,(e),-2,1,
    (-1/2e),5/2,1,(-1/2e),-5/2,1,(1/2e),5/2,1,(1/2e),-5/2,1,
    (3/2e),3/2,1,(-3/2e),-3/2,1,(3/2e),-3/2,1,(-3/2e),3/2,1,
    0,3,1,0,-3,1,(-1/4e),-1/4,1,(1/4e),1/4,1,
    (1/4e),-1/4,1,(-1/4e),1/4,1,0,-1/2,1,0,1/2,1,
    (-e),-1,1,(e),1,1,(-e),1,1,(e),-1,1,
    0,2,1,0,-2,1,(-1/2e),-1/2,1,(1/2e),1/2,1,
    (-1/2e),1/2,1,(1/2e),-1/2,1,0,1,1,0,-1,1,
    (3/2e),0,1,(-3/2e),0,1,(3/4e),9/4,1,(-3/4e),-9/4,1,
    (3/4e),-9/4,1,(-3/4e),9/4,1,(3e),0,1,(-3e),0,1,
    (-3/2e),-9/2,1,(3/2e),9/2,1,(-3/2e),9/2,1,(3/2e),-9/2,1,
    (-1/3e),0,1,(1/3e),0,1,(-1/6e),-1/2,1,(1/6e),1/2,1,
    (1/6e),-1/2,1,(-1/6e),1/2,1,(2e),0,1,(-2e),0,1,
    (-e),-3,1,(e),3,1,(-e),3,1,(e),-3,1,
    (-1/2e),0,1,(1/2e),0,1,(1/4e),3/4,1,(-1/4e),-3/4,1,
    (-1/4e),3/4,1,(1/4e),-3/4,1,(e),0,1,(-e),0,1,
    (-1/2e),-3/2,1,(1/2e),3/2,1,(-1/2e),3/2,1,(1/2e),-3/2,1,
    0,0,1,(e),1,0,(-e),1,0,0,1,0,-1,0,0,1,(e),0,-1,(e),0;

    "generating ideals I^(3) and I^2...";
    ideal I=1; ideal I3=1;
    for(int i=1;i<=(size(L) div 3);i++){
    ΨI=intersect(I,PtsIdeal(L[3*i-2],L[3*i-1],L[3*i]));
    ΨI3=intersect(I3,(PtsIdeal(L[3*i-2],L[3*i-1],L[3*i]))^3);
    Ψif((i mod 10) == 0){ string(i)+" points of "
            +string(size(L) div 3)+" in total used";}
    }
    I=std(I^2); I3=std(I3);

    "number of generators of I^(3) not in I^2: "+string(size(NF(I3,I)));

## References

  * [1] Artés, J. C., Grünbaum, B., Llibre, J.: On the number of invariant straight lines for polynomial differential systems. Pacific J. Math. 184: 207 – 230 (1998).
  * [2] Bisui, S., Grifo, E., Hà, H. T., Nguyên, T. T.: Demailly’s conjecture and the containment problem. J. Pure Appl. Algebra 226(4): Paper No. 106863, 21 (2022).
  * [3] Bocci, C., Harbourne., B.: Comparing powers and symbolic powers of ideals. J. Algebraic Geom. 19(3): 399 – 417 (2010).
  * [4] Czapliński, A., Główka, A., Malara, G., Lampa-Baczyńska, M., Łuszcz-Świdecka, P., Pokora, P., Szpond,J.: A counterexample to the containment $I^{(3)}\subset`I^{2}$ over the reals. Adv. Geom. 16(1): 77 – 82 (2016).
  * [5] Chudnovsky, G. V.: Singular points on complex hypersurfaces and multidimensional Schwarz lemma. Seminar on Number Theory, Paris 1979–80, vol. 12 Progr. Math., Birkhäuser, Boston, Mass., 29–69 (1981).
  * [6] Cook II, D., Harbourne, B., Migliore, J., Nagel, U.: Line arrangements and configurations of points with an unusual geometric property. Compos. Math. 154(10): 2150 – 2194 (2018).
  * [7] Cooper, S. M., Embree, R. J. D., Hà, H. T., Hoefel, A. H.: Symbolic powers of monomial ideals. Proc. Edinb. Math. Soc. 60(1): 39 – 55 (2017).
  * [8] Cuntz, M.: Simplicial arrangements with up to $27$ lines. Discrete Comput. Geom. 48: 682 – 701 (2012).
  * [9] Cuntz, M., Hoge, T.: Free but not recursively free arrangements. Proceedings of the American Mathematical Society 143(1): 35 – 40 (2015).
  * [10] Demailly, J.-P.: Formules de Jensen en plusieurs variables et applications arithmétiques. Bull. Soc. Math. France 110(1): 75 – 102 (1982).
  * [11] Drabkin, B., Seceleanu, A.: Singular loci of reflection arrangements and the containment problem. Math. Z. 299(1-2): 867 – 895 (2021).
  * [12] Di Marca, M., Malara, G., Oneto, A.: Unexpected curves arising from special line arrangements. J. Algebr. Comb. 51: 171 – 194 (2020).
  * [13] Dimca A.: Hyperplane Arrangements: An Introduction. Universitext, Springer, 2017.
  * [14] Dumnicki, M.: Containments of symbolic powers of ideals of generic points in $\mathbb{P}^{3}$. Proc. Amer. Math. Soc. 143(2): 513 – 530 (2015).
  * [15] Dumnicki, M., Szemberg, T., Tutaj-Gasińska, H.: Counterexamples to the $I^{(3)}\subset I^{2}$ containment. J. Algebra 393: 24 – 29 (2013).
  * [16] Dumnicki, M., Tutaj-Gasińska, H.: A containment result in $\mathbb{P}^{n}$ and the Chudnovsky conjecture. Proc. Amer. Math. Soc. 145(9): 3689 – 3694 (2017).
  * [17] Eisenbud, D., Mazur, B.: Evolutions, symbolic squares, and Fitting ideals. J. Reine Angew. Math. 488: 189 – 202 (1997).
  * [18] Ein, L., Lazarsfeld, R., Smith, K. E.: Uniform bounds and symbolic powers on smooth varieties. Invent. Math. 144(2): 241 – 252 (2001).
  * [19] Farnik, Ł., Kabat, J., Lampa-Baczyńska, M., Tutaj-Gasińska, H.: Containment problem and combinatorics. J. Algebraic Combin. 50(1): 39 – 47 (2019).
  * [20] Fouli, L., Mantero, P., Xie, Y.: Chudnovsky’s conjecture for very general points in $\mathbb{P}^{N}_{k}$. J. Algebra, 498: 211 – 227 (2018).
  * [21] Grothendieck, A. Sur La Classification Des Fibres Holomorphes Sur La Sphere De Riemann. Amer. J. Math. 79(1): 121 – 138 (1957).
  * [22] Grünbaum, B.: A catalogue of simplicial arrangements in the real projective plane. Ars Math. Contemp. 2: 1–25 (2009).
  * [23] Harbourne, B., Huneke, C.: Are symbolic powers highly evolved? J. Ramanujan Math. Soc. 28A: 247 – 266 (2013).
  * [24] Hochster, M.: Properties of Noetherian rings stable under general grade reduction. Arch. Math. (Basel) 24: 393 – 396 (1973).
  * [25] Janasz, M., Lampa-Baczyńska, M., Malara, G.: New phenomena in the containment problem for simplicial arrangements. J. Commut. Algebra 14: 571 – 581 (2022).
  * [26] Malara, G., Szpond, J.: The containment problem and a rational simplicial arrangement. Electron. Res. Announc. Math. Sci. 24: 123 -– 128 (2017).
  * [27] Malara, G., Szemberg, T., Szpond,J.: On a conjecture of Demailly and new bounds on Waldschmidt constants in $\mathbb{P}^{N}$. J. Number Theory 189: 211 – 219 (2018).
  * [28] Orlik, P., Terao, H.: Arrangements of hyperplanes, Grundlehren der math. Wiss., 300, Springer-Verlag, Berlin- Heidelberg-New York, 1992.
  * [29] Skoda, H.: Estimations $L^{2}$ pour l’opérateur $\overline{\partial}$ et applications arithmétiques. Journées sur les Fonctions Analytiques (Toulouse, 1976), Lecture Notes in Math., Vol. 578, 314–323 (1977).
  * [30] Szemberg, T., Szpond, J.: On the containment problem. Rend. Circ. Mat. Palermo (2) 66: 233 – 245 (2017).
  * [31] Terao, H.: Arrangements of hyperplanes and their freeness. I, J. Fac. Sci. Univ. Tokyo Sect. IA Math. 27(2): 293 – 312 (1980)
  * [32] Terao, H.: Generalized Exponents of a Free Arrangement of Hyperplanes and Shepherd-Todd- Brieskorn Formula. Invent. Math. 63: 159 – 179 (1981).
  * [33] Waldschmidt, M.: Propriétés arithmétiques de fonctions de plusieurs variables. II. Séminaire Pierre Lelong (Analyse) année 1975/76, Lecture Notes in Math., Vol. 578, 108–135 (1977).
  * [34] Ziegler, G.: Algebraic combinatorics of hyperplane arrangements. PhD Thesis, M.I.T., Cambridge, MA (1987).

Marek Janasz, Department of Mathematics, Pedagogical University of Krakow,
Podchora̧żych 2, PL-30-084 Kraków, Poland.
E-mail address<EMAIL_ADDRESS>