# Nonsmooth Quasistatic Modeling of Hydraulic Actuators

Ryo Kikuuwe11footnotemark: 1 44footnotemark: 4, Tomofumi Okada22footnotemark:
2 33footnotemark: 3, Hideo Yoshihara22footnotemark: 2,
Takayuki Doi33footnotemark: 3, Takao Nanjo33footnotemark: 3 and Koji
Yamashita33footnotemark: 3

###### Abstract

This article presents a quasistatic model of a hydraulic actuator driven by a
four-valve independent metering circuit. The presented model describes the
quasistatic balance between the velocity and force and that between the
flowrate and the pressure. In such balanced states, the pressure difference
across each valve determines the oil flowrate through the valve, the oil
flowrate into the actuator determines the velocity of the actuator, and the
pressures in the actuator chambers are algebraically related to the external
force acting on the actuator. Based on these relations, we derive a set of
quasistatic representations, which analytically relates the control valve
openings, the actuator velocity, and the external force. This analytical
expression is written with nonsmooth functions, of which the return values are
set-valued instead of single-valued. We also show a method of incorporating
the obtained nonsmooth quasistatic model into multibody simulators, in which a
virtual viscoelastic element is used to mimic transient responses. In
addition, the proposed model is extended to include a regeneration pipeline
and to deal with a collection of actuators driven by a single pump.

00footnotetext: Graduate School of Advanced Science and Engineering, Hiroshima
University, 1-4-1 Kagamiyama, Higashi-Hiroshima, Hiroshima 739-8527, Japan.
e-mail<EMAIL_ADDRESS>Kobelco Construction Machinery Dream-
Driven Co-Creation Research Center, Hiroshima University, 1-4-1 Kagamiyama,
Higashi-Hiroshima, Hiroshima 739-8527, Japan. 00footnotetext: Kobelco
Construction Machinery Co., Ltd., 2-2-1 Itsukaichiko, Saeki-ku, Hiroshima
731-5161, Japan.00footnotetext: Corresponding Author

> Keywords: nonsmooth hydraulics, differential-algebraic relaxation, hydraulic
> cylinders

##  1 Introduction

Control technology for construction machines requires continuing research and
development for future applications, such as remote and semi-automatic
operation. The productivity of research and development heavily depends on
simulation techniques. In particular, hydraulic actuators and hydraulic
circuits are important components of construction machines, of which the
physical behaviors need to be appropriately modeled in simulators. The
behavior of a hydraulic actuator is highly involving, depending on the oil
supply from the pump, the external forces acting on the actuator, and the
states and the characteristics of many valves in the circuit.

Many of the previous studies on the modeling of hydraulic systems assume that
the pressure is governed by first-order dynamics. The circuit is often divided
into several oil volumes, such as those in the circuit pipelines and the
actuator chambers. In each volume, the rate-of-change of the pressure is
determined by the oil flowrates in and out of the volume. Such an approach,
sometimes referred to as a lumped fluid approach, has been employed for the
controller design [32, 7, 23, 8] and for simulation purposes [34, 24].
Coupling of hydraulic systems and multibody systems have also been studied
[34, 22, 21, 20].

In the conventional model of the pressure dynamics, the rate-of-change of the
pressure is proportional to the bulk modulus divided by the volume of the oil.
The bulk modulus of the oil is usually high, and the oil volume may be small
when one needs to deal with a small segment in the pipes and when the piston
approaches either end of the cylinder. Therefore, the differential equations
representing the pressure dynamics can become numerically stiff, demanding a
small timestep size and a high computational cost for use in simulation. Some
researchers [28, 21, 13] applied the singular perturbation theory to avoid the
numerical stiffness of the governing differential equations of the pressure.
In the singular perturbation approach, the pressure dynamics is assumed to be
so fast that the steady-state pressures are quickly achieved. Following this
notion, the pressure rate-of-change in the governing differential equation is
replaced by zero, and the pressure dynamics is converted into an algebraic
constraint determining the steady-state pressure. Kiani Oshtorjani et al. [13]
have explored the applicability of this scheme to distinguish which pressure
rates-of-change can or cannot be zeroed based on the exponential stability of
the original differential equation of the pressure dynamics.

Figure 1: Hydraulic actuator and its circuit.

This article proposes a computationally-efficient modeling scheme for
hydraulic actuators, particularly focusing on the circuit structure shown in
Fig. 1. This approach focuses on the quasistatic111We use the term
‘quasistatic’ because the model does not involve the pressure dynamics but
involves the motion of the actuator and the oil.  balances established at the
steady state, at which the pressure difference across each valve determines
the oil flowrate through the valve, the oil flowrate into the actuator
determines the actuator’s velocity, the pressures in the actuator chambers
determine the force generated by the actuator, and the generated force equals
the external force acting on the actuator. This article elaborately derives an
analytical expression of the quasistatic relation between the actuator
velocity and the external force from the quasistatic representations of all
valves in the circuit. One important feature of the presented expression is
that it is nonsmooth, allowing for the set-valuedness of the pressure
especially when the valves are closed or the actuator is stopped.

The presented approach may be viewed as a full application of the singular
perturbation theory in the sense that it provides an approximate solution by
neglecting the first-order pressure dynamics. In this approach, the pressures
and flowrates are constrained algebraically, instead of through the dynamics.
In contrast to previous works [28, 21, 13], the presented nonsmooth model
allows for multiple steady states by involving the set-valuedness, and
describes the quasistatic balance between the chamber pressures and the
external force, to which the system should converge in the steady state.

In addition to the quasistatic model of Fig. 1, this article presents a method
to use the model in multibody dynamics simulation. In simulations, the
actuator model is connected with other mechanical components through a virtual
viscoelastic element to deal with the set-valuedness of the actuator model.
This approach is what Kikuuwe [15] has been referring to as a differential-
algebraic relaxation, which has been applied to many nonsmooth problems by
Kikuuwe and his colleagues [17, 18, 14, 30, 31, 16]. Some extensions toward
more complicated circuit structures are also presented.

This article is organized as follows. Section 2 shows mathematical
preliminaries, including definitions of relevant functions and some theorems
and propositions. Section 3 constructs a nonsmooth quasistatic model of the
hydraulic circuit of Fig. 1, which describes the algebraic relation among the
valve openings, the rod velocity, and the external force. Section 4 presents
an approach to incorporate the quasistatic model into multibody simulators.
Section 5 presents an extended model including a regeneration pipeline, and
Section 6 presents a model including multiple actuators driven by a single
pump. Section 7 provides some concluding remarks.

##  2 Mathematical Preliminary

###  2.1 Some Nonsmooth Functions

In this article, $\mathbb{R}$ denotes the set of all real numbers. This
article extensively uses mathematical notations of the nonsmooth system
theory, which involves set-valued functions. We use the following set-valued
functions:

$\displaystyle{\mathcal{N}}_{[A,B]}(x)$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle\left\\{\begin{array}[]{ll}[0,\infty)&\mbox{if}\ x=B\\\
0&\mbox{if}\ x\in(A,B)\\\ (-\infty,0]&\mbox{if}\ x=A\\\
\emptyset&\mbox{otherwise}\end{array}\right.$ (5)
$\displaystyle\mathrm{gsgn}(a,x,b)$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle\left\\{\begin{array}[]{ll}b&\mbox{if}\ x>0\\\
\mathrm{cl}\\{a,b\\}&\mbox{if}\ x=0\\\ a&\mbox{if}\ x<0.\end{array}\right.$
(9)

The definition (5) assumes $A<B$. Here, $\mathrm{cl}\\{a,b\\}$ stands for the
convex closure of the set $\\{a,b\\}$, being the closed set $[a,b]$ if $a\leq
b$ and $[b,a]$ if $b\leq a$. The function ${\mathcal{N}}_{{\mathcal{F}}}(x)$
is referred to as the normal cone of the set ${\mathcal{F}}$ at the point $x$.
The function $\mathrm{gsgn}$ can be seen as a generalized version of the set-
valued signum function. With the normal cone ${\mathcal{N}}$, the following
relation holds true:

$\displaystyle 0\leq x\perp y\geq 0\ \iff\ x\in-{\mathcal{N}}_{[0,\infty)}(y)\
\iff\ y\in-{\mathcal{N}}_{[0,\infty)}(x).$ (10)

Each of the above three expressions means that $x$ and $y$ are non-positive
scalars at least one of which is zero. This relation is convenient to describe
the flowrate-pressure relation at check valves. In addition, with $A<B$, the
following relation holds true:

$\displaystyle y\in{\mathcal{N}}_{[A,B]}(x)\ \iff\ x\in\mathrm{gsgn}(A,y,B).$
(11)

The following function represents the projection onto a closed set:

$\displaystyle\mathrm{proj}_{[A,B]}(x)$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle\max(A,\min(B,x))$ (12)

where $A<B$. The normal cone and the projection have the following relation:

$\displaystyle a-x\in{\mathcal{N}}_{{\mathcal{A}}}(x)\ \iff\
x=\mathrm{proj}_{{\mathcal{A}}}(a),$ (13)

which has been shown in previuos publications [5, Proposition 2][2, Section
A.3]. This article uses the following theorem:

######  Theorem 1.

Let $x\in\mathbb{R}$ and let ${\mathcal{A}}\subset\mathbb{R}$ be a closed
subset of $\mathbb{R}$. Let $f:\mathbb{R}\to\mathbb{R}$ be a strictly
decreasing function. Let $x_{f}\in\mathbb{R}$ satisfy $f(x_{f})=0$. Then, the
following statement holds true:

$\displaystyle f(x)\in{\mathcal{N}}_{{\mathcal{A}}}(x)\ \iff\
x=\mathrm{proj}_{{\mathcal{A}}}(x_{f}).$ (14)

###### Proof.

Because $f$ is strictly decreasing, $x<x_{f}\iff f(x)<0$, $x>x_{f}\iff
f(x)>0$, and $x=x_{f}\iff f(x)=0$ are satisfied. This means that
$f(x)\in{\mathcal{N}}_{{\mathcal{A}}}(x)\iff
x_{f}-x\in{\mathcal{N}}_{{\mathcal{A}}}(x)$ and (13) implies that it is
equivalent to $x=\mathrm{proj}_{{\mathcal{A}}}(x_{f})$. ∎

###  2.2 Some Smooth Functions

The following single-valued functions are used in the article:

$\displaystyle{\mathcal{S}}(x)$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle\mathrm{sgn}(x)x^{2}$ (15) $\displaystyle{\mathcal{R}}(x)$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle\mathrm{sgn}(x)\sqrt{|x|}$ (16) $\displaystyle\psi(u_{1},u_{2})$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle\left\\{\begin{array}[]{ll}0&\mbox{if }\ u_{1}=u_{2}=0\\\
\dfrac{u_{1}^{2}u_{2}^{2}}{u_{1}^{2}+u_{2}^{2}}&\mbox{otherwise.}\end{array}\right.$
(19)

The functions ${\mathcal{R}}$ and ${\mathcal{S}}$ are strictly increasing
continuous functions satisfying
${\mathcal{R}}({\mathcal{S}}(x))={\mathcal{S}}({\mathcal{R}}(x))=x$. The
function $\psi$ is differentiable everywhere.

Section 4 will use the following functions:

$\displaystyle\Phi_{A}(b,c,a)$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle-\mathrm{sgn}(c)\dfrac{\sqrt{a^{2}b^{2}+4a|c|}-ab}{2}$ (20)
$\displaystyle\Phi_{B}(b,c,a_{0},a_{1},x_{1})$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle\left\\{\begin{array}[]{lr}-\dfrac{\sqrt{a_{0}^{2}(a_{1}b+2x_{1})^{2}+4a_{0}(a_{0}+a_{1})(a_{1}c-x_{1}^{2})}-a_{0}(a_{1}b+2x_{1})}{2(a_{0}+a_{1})}&\\\\[12.0pt]
&\mbox{if }\left(-a_{0}(bx+c)\leq{\mathcal{S}}(x_{1})\leq
0\right)\,\vee\,\left(0\leq{\mathcal{S}}(x_{1})\leq a_{1}c\right)\\\\[12.0pt]
\dfrac{\sqrt{a_{0}^{2}(a_{1}b-2x_{1})^{2}-4a_{0}(a_{0}+a_{1})(a_{1}c+x_{1}^{2})}-a_{0}(a_{1}b-2x_{1})}{2(a_{0}+a_{1})}&\\\\[12.0pt]
&\mbox{if }\left(0\leq{\mathcal{S}}(x_{1})\leq-
a_{0}(bx+c)\right)\,\vee\,\left(a_{1}c\leq{\mathcal{S}}(x_{1})\leq
0\right)\\\\[12.0pt]
-\dfrac{2\sqrt{a_{0}}(a_{1}c+x_{1}^{2})}{\sqrt{a_{0}}(a_{1}b-2x_{1})+\sqrt{a_{0}(a_{1}b-2x_{1})^{2}+4(a_{1}-a_{0})(a_{1}c+x_{1}^{2})}}&\\\\[12.0pt]
&\mbox{if }{\mathcal{S}}(x_{1})\leq\min(0,-a_{0}(bx+c),a_{1}c)\\\\[12.0pt]
\dfrac{2\sqrt{a_{0}}(x_{1}^{2}-a_{1}c)}{\sqrt{a_{0}}(a_{1}b+2x_{1})+\sqrt{a_{0}(a_{1}b+2x_{1})^{2}+4(a_{1}-a_{0})(x_{1}^{2}-a_{1}c)}}&\\\\[12.0pt]
&\mbox{if
}\max(0,-a_{0}(bx_{1}+c),a_{1}c)\leq{\mathcal{S}}(x_{1}).\end{array}\right.$
(29)

With these functions, we have the following propositions:

######  Proposition 1.

Let $a\geq 0$ and $b>0$. Then, the following statement holds true:

$\displaystyle x=\Phi_{A}(b,c,a)\ \iff\ \lim_{\tilde{a}\to
a}\dfrac{{\mathcal{S}}(x)}{\tilde{a}}+bx+c=0.$ (30)

######  Proposition 2.

Let $a_{0}\geq 0$, $a_{1}\geq 0$, and $b>0$, and let $a_{0}^{2}+a_{1}^{2}>0$
be satisfied. Then, the following statement holds true:

$\displaystyle x=\Phi_{B}(b,c,a_{0},a_{1},x_{1})\ \iff\ \lim_{\tilde{a}_{0}\to
a_{0}}\dfrac{{\mathcal{S}}(x)}{\tilde{a}_{0}}+\lim_{\tilde{a}_{1}\to
a_{1}}\dfrac{{\mathcal{S}}(x-x_{1})}{\tilde{a}_{1}}+bx+c=0.$ (31)

The proofs of these propositions can be obtained through tedious but
straightforward derivations.

We also use the following theorems, of which the proofs are rather trivial:

######  Theorem 2.

Let $f:\mathbb{R}\to\mathbb{R}$ and $g:\mathbb{R}\to\mathbb{R}$ be strictly
decreasing functions and let $g(0)\leq f(0)$ be satisfied. Let $x_{f}$ and
$x_{g}$ satisfy $b=f(x_{f})$ and $b=g(x_{g})$ where $b\in\mathbb{R}$. Then,
the following statement holds true:

$\displaystyle b\in\mathrm{gsgn}(f(x),x,g(x))\ \iff\
x=\left\\{\begin{array}[]{ll}x_{f}&\mbox{if }b>f(0)\\\ x_{g}&\mbox{if
}b<g(0)\\\ 0&\mbox{otherwise}.\\\ \end{array}\right.$ (35)

######  Theorem 3.

Let $f_{1}:\mathbb{R}\to\mathbb{R}$ and $f_{2}:\mathbb{R}\to\mathbb{R}$ be
strictly decreasing functions. Let $x_{i}$ satisfy $f_{i}(x)=0$ for
$i\in\\{1,2\\}$. Then, the following two statements hold true:

$\displaystyle\min(f_{1}(x),f_{2}(x))=0$ $\displaystyle\iff$ $\displaystyle
x=\min(x_{1},x_{2})$ (36) $\displaystyle\max(f_{1}(x),f_{2}(x))=0$
$\displaystyle\iff$ $\displaystyle x=\max(x_{1},x_{2}).$ (37)

##  3 Nonsmooth Quasistatic Model

This section considers the hydraulic circuit illustrated in Fig. 1, which is a
four-valve independent metering circuit to drive a double-acting hydraulic
actuator. The actuator has two chambers separated by the piston, and the
motion of the piston is extracted as the motion of the rod, which applies
forces to external objects. We are interested in the quasistatic relation
among the rod velocity $v$, positive when the rod is extending, the external
force $f$, positive when it is compressing the rod, and the opening ratios of
the valves. This hydraulic circuit is similar to those studied in, e.g., [27,
25, 26, 9, 6], where the quasistatic relations are also considered. Our main
contribution lies in an elaborate analytical representation of the whole
circuit, which is rather complicated than those in previous studies, using the
nonsmooth formalism to deal with relief valves and check valves.

This article uses the terminology for linear hydraulic actuators (i.e.,
hydraulic cylinders), but the presented approach is applicable also to rotary
hydraulic actuators by replacing the velocity and the external force by the
angular velocity and the torque, respectively.

###  3.1 Quasistatic Relations

In Fig. 1, $Q_{*}$ and $P_{*}$ denote the flowrates and the pressure at each
point. The pump provides the flowrate $Q$ to the circuit via a pump check
valve and two of the four main control valves. These control valves are
connected to the head-side and the rod-side chambers of the actuator. The
chambers are also connected to the tank with the zero pressure via the other
two control valves. Each chamber of the actuator also connects to the tank
through a parallel combination of a relief valve and a check valve, which are
named as indicated in the figure. There is another control valve, referred to
as a bleed valve, and it leads to the tank in parallel to a relief valve,
referred to as a pump relief valve.

Figure 2: Flowrate-pressure relations at: (a) one of main control valves
(39h), (b) the bleed and the pump relief valves (39l), (c) the pump check
valve (39m), and (d) the head-side relief and the suction check valves (39n).

The degrees of the opening of the valves are represented by dimensionless
variables $u_{*}\in[0,1]$ ($*\in\\{ph,tr,pr,th,b\\}$), which are the ratios of
the valve opening areas to their maximum values. The control valves are
manipulated by a controller that accepts external commands, which are given
through, e.g., operation levers of an excavator. When the external command is
to stop the actuator, all the four main valves are closed. When the external
command is to move the actuator in the positive direction (i.e., extend the
rod), both or either of $u_{ph}$ and $u_{tr}$ are set positive and $u_{th}$
and $u_{pr}$ are set zero. When the operator’s command is to move the actuator
in the negative direction (i.e, retract the rod), $u_{ph}$ and $u_{tr}$ are
set zero and both or either of $u_{th}$ and $u_{pr}$ are set positive. Based
on this idea, we assume that the vector
${\boldsymbol{u}}\stackrel{{\scriptstyle\Delta}}{{=}}[u_{ph},u_{tr},u_{pr},u_{th},u_{b}]^{T}$
always belongs to either of the following three subsets:

$\displaystyle{\boldsymbol{\mathcal{U}}}_{0}$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle\\{{\boldsymbol{u}}\in\mathbb{R}^{5}\,|\,u_{1}=u_{2}=u_{3}=u_{4}=0\,\wedge\,u_{5}>0\\}$
(38a) $\displaystyle{\boldsymbol{\mathcal{U}}}_{+}$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle\\{{\boldsymbol{u}}\in\mathbb{R}^{5}\,|\,u_{1}^{2}+u_{2}^{2}>0\,\wedge\,u_{3}=u_{4}=0\,\wedge\,u_{5}\geq
0\\}$ (38b) $\displaystyle{\boldsymbol{\mathcal{U}}}_{-}$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle\\{{\boldsymbol{u}}\in\mathbb{R}^{5}\,|\,u_{1}=u_{2}=0\,\wedge\,u_{3}^{2}+u_{4}^{2}>0\,\wedge\,u_{5}\geq
0\\}$ (38c)

where $u_{i}$ ($i\in\\{1,\cdots,5\\}$) stands for the $i$th element of the
vector ${\boldsymbol{u}}$.

Now, let us make an exhaustive list of algebraic relations among the
pressures, the flowrates, the external force, and the actuator velocity at the
steady state. First, according to the principle of mass conservation, one can
see that the following relations hold true at junctions in the circuit:

$\displaystyle Q_{hn}$ $\displaystyle=$ $\displaystyle Q_{ph}-Q_{th}$ (39a)
$\displaystyle Q_{rn}$ $\displaystyle=$ $\displaystyle Q_{pr}-Q_{tr}$ (39b)
$\displaystyle Q_{h}$ $\displaystyle=$ $\displaystyle Q_{hn}-Q_{ho}$ (39c)
$\displaystyle Q_{r}$ $\displaystyle=$ $\displaystyle Q_{rn}-Q_{ro}$ (39d)
$\displaystyle Q_{t}$ $\displaystyle=$ $\displaystyle Q_{th}+Q_{tr}$ (39e)
$\displaystyle Q_{p}$ $\displaystyle=$ $\displaystyle Q_{ph}+Q_{pr}$ (39f)
$\displaystyle Q_{b}$ $\displaystyle=$ $\displaystyle Q-Q_{p}.$ (39g)

Second, let us focus on the control valves. As indicated in Fig. 1, $P_{h}$
and $P_{r}$ are the internal pressures of the head- and rod-side chambers,
respectively, $P_{c}$ is the pressure at the check valve connected to the
pump, and $P$ is the pressure at the outlet of the pump. According to the
conventional orifice model [4, 7], we can assume that the following flowrate-
pressure relations are satisfied:

$\displaystyle Q_{ph}$ $\displaystyle=$ $\displaystyle
c_{ph}u_{ph}{\mathcal{R}}(P_{c}-P_{h})$ (39h) $\displaystyle Q_{th}$
$\displaystyle=$ $\displaystyle c_{th}u_{th}{\mathcal{R}}(P_{h})$ (39i)
$\displaystyle Q_{pr}$ $\displaystyle=$ $\displaystyle
c_{pr}u_{pr}{\mathcal{R}}(P_{c}-P_{r})$ (39j) $\displaystyle Q_{tr}$
$\displaystyle=$ $\displaystyle c_{tr}u_{tr}{\mathcal{R}}(P_{r})$ (39k)
$\displaystyle Q_{b}$ $\displaystyle\in$ $\displaystyle
c_{b}u_{b}{\mathcal{R}}(P)+{\mathcal{N}}_{(-\infty,P_{M}]}(P).$ (39l)

Fig. 2(a) illustrates the relation (39h). The coefficients $c_{*}$ are defined
as $c_{*}=C_{*}a_{*}\sqrt{2/\rho}$ ($*\in\\{ph,tr,pr,th,b\\}$) where $C_{*}$
is the dimensionless coefficient named a discharge coefficient [19], which is
typically around $0.6$ or $0.7$ [33, 29], $a_{*}$ is the maximum opening area
(m2) of the valve, and $\rho$ is the mass density (kg/m3) of the oil. Equation
(39l) represents the combined effect of the relief valve and the bleed valve,
which limits the pump pressure up to $P_{M}$ as illustrated in Fig. 2(b).

Third, let us consider the check valves and the relief valves. The check valve
connected to the pump imposes the following constraint:

$\displaystyle Q_{p}$ $\displaystyle\in$
$\displaystyle-{\mathcal{N}}_{[0,\infty)}(P_{c}-P),$ (39m)

which means that the flowrate $Q_{p}$ is zero when $P_{c}-P>0$, as illustrated
in Fig. 2(c). The effects of the relief valves and the suction check valves
connected to the chambers are written as follows:

$\displaystyle Q_{ho}$ $\displaystyle\in$
$\displaystyle{\mathcal{N}}_{[0,P_{hM}]}(P_{h})$ (39n) $\displaystyle Q_{ro}$
$\displaystyle\in$ $\displaystyle{\mathcal{N}}_{[0,P_{rM}]}(P_{r}).$ (39o)

Here, $P_{hM}$ and $P_{rM}$ are the pressure limits of the relief valves
connected to the head- and rod-side chambers, respectively. These expressions
mean that $P_{h}$ and $P_{r}$ are always in the ranges of $[0,P_{hM}]$ and
$[0,P_{rM}]$, respectively, as illustrated in Fig. 2(d). When the pressure
reaches the upper limit, the oil flows into the tank. When the pressure
reaches zero, the oil is drawn into the chamber from the tank.

Lastly, we discuss the actuator. Let $A_{h}$ and $A_{r}$ be the cross-
sectional areas of the head- and rod-side chambers, respectively. Then, at the
steady state where the rod inertia can be neglected, the constraints imposed
by the actuator can be written as follows:

$\displaystyle v$ $\displaystyle=$ $\displaystyle Q_{h}/A_{h}$ (39p)
$\displaystyle v$ $\displaystyle=$ $\displaystyle-Q_{r}/A_{r}$ (39q)
$\displaystyle f$ $\displaystyle=$ $\displaystyle A_{h}P_{h}-A_{r}P_{r}.$
(39r)

If one deals with rotary actuators, both $A_{r}$ and $A_{h}$ (measured in m2)
should be replaced by the volume displacement per one radian of rotation,
which is measured in m3/rad.

In conclusion, now we have 18 algebraic constraints in (39) and 19 variables,
which are listed as follows

  * •

$\\{Q_{ho},Q_{ro},Q_{hn},Q_{rn},Q_{h},Q_{r},Q_{ph},Q_{pr},Q_{th},Q_{tr},Q_{t},Q_{p},Q_{b}\\}$:
13 flowrate values

  * •

$\\{P,P_{h},P_{r},P_{c}\\}$: 4 pressure values

  * •

$\\{f,v\\}$: the external force to the rod and the velocity of the rod.

######  Remark 1.

The presented formalism can be said to be close to the classical hydraulic-
electric analogy [12, 10], which replaces the pressure and the flowrate by the
voltage and the current, respectively. Considering that the presented approach
involves the nonsmoothness, it is also analogous to the nonsmooth electronics
[3, 11, 1]. In this analogy, a check valve, for example, corresponds to an
ideal diode [1].

###  3.2 Normalized Representations

For the convenience of derivation, we now normalize some quantities in the
following manner:

$\displaystyle q_{*}\stackrel{{\scriptstyle\Delta}}{{=}}Q_{*}/A_{h}\
(*\in\\{ho,hn,h,ph,th\\}),\quad\hat{u}_{*}\stackrel{{\scriptstyle\Delta}}{{=}}c_{*}u_{*}/A_{h}^{3/2}\
(*\in\\{ph,th\\})$ (40) $\displaystyle
q_{*}\stackrel{{\scriptstyle\Delta}}{{=}}Q_{*}/A_{r}\
(*\in\\{ro,rn,r,pr,tr\\}),\quad\hat{u}_{*}\stackrel{{\scriptstyle\Delta}}{{=}}c_{*}u_{*}/A_{r}^{3/2}\
(*\in\\{pr,tr\\})$ (41) $\displaystyle
F_{*}\stackrel{{\scriptstyle\Delta}}{{=}}P_{*}A_{*},\quad F_{*M}=P_{*M}A_{*}\
(*\in\\{h,r\\})$ (42) $\displaystyle
U_{b}\stackrel{{\scriptstyle\Delta}}{{=}}c_{b}u_{b}.$ (43)

The regularized input vector is defined as
$\hat{{\boldsymbol{u}}}=[\hat{u}_{ph},\hat{u}_{tr},\hat{u}_{pr},\hat{u}_{th},U_{b}]^{T}$.
By using these definitions, (39) can be rewritten as follows:

$\displaystyle q_{hn}=q_{ph}-q_{th}$ (44a) $\displaystyle
q_{rn}=q_{pr}-q_{tr}$ (44b) $\displaystyle q_{h}=q_{hn}-q_{ho}$ (44c)
$\displaystyle q_{r}=q_{rn}-q_{ro}$ (44d) $\displaystyle
Q_{t}=A_{h}q_{th}+A_{r}q_{tr}$ (44e) $\displaystyle
Q_{p}=A_{h}q_{ph}+A_{r}q_{pr}$ (44f) $\displaystyle Q_{b}=Q-Q_{p}$ (44g)
$\displaystyle q_{ph}=\hat{u}_{ph}{\mathcal{R}}(A_{h}P_{c}-F_{h})$ (44h)
$\displaystyle q_{th}=\hat{u}_{th}{\mathcal{R}}(F_{h})$ (44i) $\displaystyle
q_{pr}=\hat{u}_{pr}{\mathcal{R}}(A_{r}P_{c}-F_{r})$ (44j) $\displaystyle
q_{tr}=\hat{u}_{tr}{\mathcal{R}}(F_{r})$ (44k) $\displaystyle Q_{b}\in
U_{b}{\mathcal{R}}(P)+{\mathcal{N}}_{(-\infty,P_{M}]}(P)$ (44l) $\displaystyle
Q_{p}\in{\mathcal{N}}_{(-\infty,P_{c}]}(P)$ (44m) $\displaystyle
q_{ho}\in{\mathcal{N}}_{[0,F_{hM}]}(F_{h})$ (44n) $\displaystyle
q_{ro}\in{\mathcal{N}}_{[0,F_{rM}]}(F_{r})$ (44o) $\displaystyle v=q_{h}$
(44p) $\displaystyle v=-q_{r}$ (44q) $\displaystyle f=F_{h}-F_{r}.$ (44r)

Note that the equivalence between (39m) and (44m) can be derived from the
definition of the normal cone.

Now, let us attempt to reduce the number of variables. Substituting (44m) and
(44g) into (44l) yields

$\displaystyle Q$ $\displaystyle\in$
$\displaystyle{\mathcal{N}}_{(-\infty,P_{c}]}(P)+U_{b}{\mathcal{R}}(P)+{\mathcal{N}}_{(-\infty,P_{M}]}(P)$
(45) $\displaystyle=$ $\displaystyle
U_{b}{\mathcal{R}}(P)+{\mathcal{N}}_{(-\infty,\min(P_{c},P_{M})]}(P),$

which is equivalent to

$\displaystyle P$ $\displaystyle=$
$\displaystyle\min(P_{c},P_{M},{\mathcal{S}}(Q/U_{b}))$ (46)

because of Theorem 1. Since $q_{h}=-q_{r}=v$, we can eliminate the two
variables $q_{h}$ and $q_{r}$. Moreover, $q_{ho}$, $q_{ro}$, $q_{hn}$,
$q_{rn}$, $q_{ph}$, $q_{th}$, $q_{pr}$, $q_{tr}$, and $Q_{t}$ can also be
eliminated. This gives the following five equations for the four variables
$\\{P,P_{c},F_{h},F_{r},f\\}$:

$\displaystyle-v+\hat{u}_{ph}{\mathcal{R}}(A_{h}P_{c}-F_{h})-\hat{u}_{th}{\mathcal{R}}(F_{h})\in{\mathcal{N}}_{[0,F_{hM}]}(F_{h})$
(47a) $\displaystyle
v+\hat{u}_{pr}{\mathcal{R}}(A_{r}P_{c}-F_{r})-\hat{u}_{tr}{\mathcal{R}}(F_{r})\in{\mathcal{N}}_{[0,F_{rM}]}(F_{r})$
(47b) $\displaystyle
A_{h}\hat{u}_{ph}{\mathcal{R}}(A_{h}P_{c}-F_{h})+A_{r}\hat{u}_{pr}{\mathcal{R}}(A_{r}P_{c}-F_{r})\in{\mathcal{N}}_{(-\infty,P_{c}]}(P)$
(47c) $\displaystyle Q\in
U_{b}{\mathcal{R}}(P)+A_{h}\hat{u}_{ph}{\mathcal{R}}(A_{h}P_{c}-F_{h})+A_{r}\hat{u}_{pr}{\mathcal{R}}(A_{r}P_{c}-F_{r})+{\mathcal{N}}_{(-\infty,P_{M}]}(P)$
(47d) $\displaystyle f=F_{h}-F_{r}.$ (47e)

###  3.3 Main Result: Nonsmooth Quasistatic Map from $v$ to $f$

Now we derive the relation between $f$ and $v$ from (47). If
$\hat{{\boldsymbol{u}}}\in{\boldsymbol{\mathcal{U}}}_{0}$, (47) reduces to the
following:

$\displaystyle-v$ $\displaystyle\in$
$\displaystyle{\mathcal{N}}_{[0,F_{hM}]}(F_{h})$ (48a) $\displaystyle v$
$\displaystyle\in$ $\displaystyle{\mathcal{N}}_{[0,F_{rM}]}(F_{r})$ (48b)
$\displaystyle 0$ $\displaystyle\in$
$\displaystyle{\mathcal{N}}_{(-\infty,P_{c}]}(P)$ (48c) $\displaystyle Q$
$\displaystyle\in$ $\displaystyle
U_{b}{\mathcal{R}}(P)+{\mathcal{N}}_{(-\infty,P_{M}]}(P)$ (48d) $\displaystyle
f$ $\displaystyle=$ $\displaystyle F_{h}-F_{r},$ (48e)

from which

$\displaystyle F_{h}\in\mathrm{gsgn}(F_{hM},v,0),\quad
F_{r}\in\mathrm{gsgn}(0,v,F_{rM})$ (49)

and $f\in\mathrm{gsgn}(F_{hM},v,-F_{rM})$ can be derived by using (11). This
means that $f$ can take any values between $F_{hM}$ and $-F_{rM}$ when $v=0$,
which is consistent with the fact that, when all the main control valves are
closed, the cylinder holds its position by producing the reaction force
against the external force as long as the relief valves are closed.

If $\hat{{\boldsymbol{u}}}\in{\boldsymbol{\mathcal{U}}}_{+}$, (47) reduces to
the following:

$\displaystyle-v+\hat{u}_{ph}{\mathcal{R}}(A_{h}P_{c}-F_{h})\in{\mathcal{N}}_{[0,F_{hM}]}(F_{h})$
(50a) $\displaystyle
v-\hat{u}_{tr}{\mathcal{R}}(F_{r})\in{\mathcal{N}}_{[0,F_{rM}]}(F_{r})$ (50b)
$\displaystyle
A_{h}\hat{u}_{ph}{\mathcal{R}}(A_{h}P_{c}-F_{h})\in{\mathcal{N}}_{(-\infty,P_{c}]}(P)$
(50c) $\displaystyle Q\in
U_{b}{\mathcal{R}}(P)+A_{h}\hat{u}_{ph}{\mathcal{R}}(A_{h}P_{c}-F_{h})+{\mathcal{N}}_{(-\infty,P_{M}]}(P)$
(50d) $\displaystyle f=F_{h}-F_{r}.$ (50e)

Here, (50a), (50b) and (50c) can be rewritten as

$\displaystyle
F_{h}=\mathrm{proj}_{[0,F_{hM}]}(A_{h}P_{c}-{\mathcal{S}}(v)/\hat{u}_{ph}^{2})$
(51) $\displaystyle
F_{r}=\mathrm{proj}_{[0,F_{rM}]}({\mathcal{S}}(v)/\hat{u}_{tr}^{2})$ (52)
$\displaystyle P_{c}=\max(P,F_{h}/A_{h}),$ (53)

respectively, because of Theorem 1. Substituting (53) into (51) and (50d)
results in:

$\displaystyle F_{h}$ $\displaystyle=$
$\displaystyle\mathrm{proj}_{[0,F_{hM}]}(\max(A_{h}P,F_{h})-{\mathcal{S}}(v)/\hat{u}_{ph}^{2}).$
(54a) $\displaystyle Q$ $\displaystyle\in$ $\displaystyle
U_{b}{\mathcal{R}}(P)+A_{h}\hat{u}_{ph}{\mathcal{R}}(\max(A_{h}P-F_{h},0))+{\mathcal{N}}_{(-\infty,P_{M}]}(P),$
(54b)

respectively. Now $F_{r}$ is obtained by (52), and thus we focus on obtaining
$F_{h}$ from (54). If $v<0$, (54a) implies $F_{h}=F_{hM}$. If $v=0$, (54a)
implies that $\min(A_{h}P,F_{hM})\leq F_{h}\leq F_{hM}$. If $v=0$ and
$F_{h}<F_{hM}$, (54a) implies $F_{h}\geq A_{h}P$ and substituting it into
(54b) yields $P=\min(P_{M},Q^{2}/U_{b}^{2})$. Therefore, if $v=0$, the
condition (54) implies the following:

$\displaystyle F_{h}$ $\displaystyle\in$
$\displaystyle[\min(F_{hM},A_{h}P_{M},A_{h}Q^{2}/U_{b}^{2}),F_{hM}].$ (55)

If $v>0$, (54a) implies $A_{h}P>F_{h}$ and thus (54) can be rewritten as
follows:

$\displaystyle F_{h}$ $\displaystyle=$
$\displaystyle\mathrm{proj}_{[0,F_{hM}]}(A_{h}P-{\mathcal{S}}(v)/\hat{u}_{ph}^{2}).$
(56a) $\displaystyle Q$ $\displaystyle\in$ $\displaystyle
U_{b}{\mathcal{R}}(P)+A_{h}\hat{u}_{ph}{\mathcal{R}}(A_{h}P-\mathrm{proj}_{[0,F_{hM}]}(A_{h}P-{\mathcal{S}}(v)/\hat{u}_{ph}^{2}))+{\mathcal{N}}_{(-\infty,P_{M}]}(P).\quad$
(56b)

The relation (56b) implies that $P$ is a strictly increasing function of $Q$
that saturates at $P=P_{M}$. As long as
$A_{h}P-{\mathcal{S}}(v)/\hat{u}_{ph}^{2}\in[0,F_{hM}]$, $P$ is written as

$\displaystyle
P=\min\left(P_{M},-\dfrac{A_{h}^{2}}{U_{b}^{2}}{\mathcal{S}}\left(v-\dfrac{Q}{A_{h}}\right)\right).$
(57)

Outside this range, $P$ is expressed as another function of $Q$ that strictly
increases with respect to $Q$ and continuously connects to (57). The value of
$P$ outside the range does not influence $F_{h}$ because $F_{h}$ is saturated
at $0$ and $F_{hM}$ due to (56a). Therefore, if $v>0$, the condition (54)
implies

$\displaystyle F_{h}$ $\displaystyle=$
$\displaystyle\mathrm{proj}_{[0,F_{hM}]}\left(\min\left(A_{h}P_{M},-\dfrac{A_{h}^{3}}{U_{b}^{2}}{\mathcal{S}}\left(v-\dfrac{Q}{A_{h}}\right)\right)-\dfrac{{\mathcal{S}}(v)}{\hat{u}_{ph}^{2}}\right),$
(58)

which is obtained by substituting (57) into (56a). Unifying the three cases,
i.e., $F_{h}=F_{hM}$ for $v<0$, (55) for $v=0$, and (58) for $v>0$, we have
the following expression for the case
$\hat{{\boldsymbol{u}}}\in{\boldsymbol{\mathcal{U}}}_{+}$:

$\displaystyle F_{h}$ $\displaystyle\in$
$\displaystyle\mathrm{gsgn}\left(F_{hM},v,\mathrm{proj}_{[0,F_{hM}]}\left(\min\left(A_{h}P_{M},-\dfrac{A_{h}^{3}}{U_{b}^{2}}{\mathcal{S}}\left(v-\dfrac{Q}{A_{h}}\right)\right)-\dfrac{{\mathcal{S}}(v)}{\hat{u}_{ph}^{2}}\right)\right).$
(59)

If $\hat{{\boldsymbol{u}}}\in{\boldsymbol{\mathcal{U}}}_{-}$, we have

$\displaystyle
F_{h}=\mathrm{proj}_{[0,F_{hM}]}(-{\mathcal{S}}(v)/\hat{u}_{th}^{2})$ (60)

for any $v$ in the same manner as (52). Combining (60) for
$\hat{{\boldsymbol{u}}}\in{\boldsymbol{\mathcal{U}}}_{-}$, (59) for
$\hat{{\boldsymbol{u}}}\in{\boldsymbol{\mathcal{U}}}_{+}$ and (49) for
$\hat{{\boldsymbol{u}}}\in{\boldsymbol{\mathcal{U}}}_{0}$, we have

$\displaystyle F_{h}$ $\displaystyle\in$
$\displaystyle\Gamma_{h}(v)\stackrel{{\scriptstyle\Delta}}{{=}}\mathrm{gsgn}\left(\Gamma_{h-}(v),v,\Gamma_{h+}(v)\right)$
(61a) where $\displaystyle\Gamma_{h+}(v)$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle\mathrm{proj}_{[0,F_{hM}]}\left(\min\left(A_{h}P_{M},-\dfrac{A_{h}^{3}}{U_{b}^{2}}{\mathcal{S}}\left(v-\dfrac{Q}{A_{h}}\right)\right)-\dfrac{{\mathcal{S}}(v)}{\hat{u}_{ph}^{2}}\right)$
(61b) $\displaystyle\Gamma_{h-}(v)$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle\mathrm{proj}_{[0,F_{hM}]}\left(-\dfrac{{\mathcal{S}}(v)}{\hat{u}_{th}^{2}}\right)$
(61c)

for any $\hat{{\boldsymbol{u}}}$. In the same manner, $F_{r}$ for any
$\hat{{\boldsymbol{u}}}$ is obtained as follows:

$\displaystyle F_{r}$ $\displaystyle\in$
$\displaystyle\Gamma_{r}(v)\stackrel{{\scriptstyle\Delta}}{{=}}\mathrm{gsgn}\left(\Gamma_{r-}(v),v,\Gamma_{r+}(v)\right)$
(62a) where $\displaystyle\Gamma_{r+}(v)$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle\mathrm{proj}_{[0,F_{rM}]}\left(\dfrac{{\mathcal{S}}(v)}{\hat{u}_{tr}^{2}}\right)$
(62b) $\displaystyle\Gamma_{r-}(v)$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle\mathrm{proj}_{[0,F_{rM}]}\left(\min\left(A_{r}P_{M},\dfrac{A_{r}^{3}}{U_{b}^{2}}{\mathcal{S}}\left(v+\dfrac{Q}{A_{r}}\right)\right)+\dfrac{{\mathcal{S}}(v)}{\hat{u}_{pr}^{2}}\right).$
(62c)

By using (61) and (62), the relation between $v$ and $f$ is obtained as
follows:

$\displaystyle f$ $\displaystyle\in$
$\displaystyle\Gamma(v)\stackrel{{\scriptstyle\Delta}}{{=}}\Gamma_{h}(v)-\Gamma_{r}(v).$
(63)

For the convenience of further derivations, we can also write the set-valued
map $\Gamma(v)$ in the following form:

$\displaystyle\Gamma(v)$ $\displaystyle=$
$\displaystyle\mathrm{gsgn}(\Gamma_{-}(v),v,\Gamma_{+}(v))$ (64a) where
$\displaystyle\Gamma_{+}(v)$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle\Gamma_{h+}(v)-\Gamma_{r+}(v)$ (64c) $\displaystyle=$
$\displaystyle\max(\min(\max(\Gamma_{+0a}(v),\Gamma_{+0b}(v)),\max(\Gamma_{+1a}(v),\Gamma_{+1b}(v)),$
$\displaystyle\quad\max(\Gamma_{+2a}(v),\Gamma_{+2b}(v))),\Gamma_{+3}(v),-F_{rM})\quad$
$\displaystyle\Gamma_{-}(v)$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle\Gamma_{h-}(v)-\Gamma_{r-}(v)$ (64e) $\displaystyle=$
$\displaystyle\min(\max(\min(\Gamma_{-0a}(v),\Gamma_{-0b}(v)),\min(\Gamma_{-1a}(v),\Gamma_{-1b}(v)),$
$\displaystyle\quad\min(\Gamma_{-2a}(v),\Gamma_{-2b}(v))),\Gamma_{-3}(v),F_{hM})\quad$
$\displaystyle\Gamma_{+0a}(v)$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$ $\displaystyle
F_{hM}-\dfrac{{\mathcal{S}}(v)}{\hat{u}_{tr}^{2}}$ (64f)
$\displaystyle\Gamma_{+0b}(v)$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$ $\displaystyle
F_{hM}-F_{rM}$ (64g) $\displaystyle\Gamma_{+1a}(v)$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle-\dfrac{A_{h}^{3}}{U_{b}^{2}}{\mathcal{S}}\left(v-\dfrac{Q}{A_{h}}\right)-\dfrac{{\mathcal{S}}(v)}{\hat{u}_{ph}^{2}}-\dfrac{{\mathcal{S}}(v)}{\hat{u}_{tr}^{2}}$
(64h) $\displaystyle\Gamma_{+1b}(v)$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle-\dfrac{A_{h}^{3}}{U_{b}^{2}}{\mathcal{S}}\left(v-\dfrac{Q}{A_{h}}\right)-\dfrac{{\mathcal{S}}(v)}{\hat{u}_{ph}^{2}}-F_{rM}$
(64i) $\displaystyle\Gamma_{+2a}(v)$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$ $\displaystyle
A_{h}P_{M}-\dfrac{{\mathcal{S}}(v)}{\hat{u}_{ph}^{2}}-\dfrac{{\mathcal{S}}(v)}{\hat{u}_{tr}^{2}}$
(64j) $\displaystyle\Gamma_{+2b}(v)$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$ $\displaystyle
A_{h}P_{M}-\dfrac{{\mathcal{S}}(v)}{\hat{u}_{ph}^{2}}-F_{rM}$ (64k)
$\displaystyle\Gamma_{+3}(v)$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle-\dfrac{{\mathcal{S}}(v)}{\hat{u}_{tr}^{2}}$ (64l)
$\displaystyle\Gamma_{-0a}(v)$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$ $\displaystyle-
F_{rM}-\dfrac{{\mathcal{S}}(v)}{\hat{u}_{th}^{2}}$ (64m)
$\displaystyle\Gamma_{-0b}(v)$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$ $\displaystyle-
F_{rM}+F_{hM}$ (64n) $\displaystyle\Gamma_{-1a}(v)$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle-\dfrac{A_{r}^{3}}{U_{b}^{2}}{\mathcal{S}}\left(v+\dfrac{Q}{A_{r}}\right)-\dfrac{{\mathcal{S}}(v)}{\hat{u}_{pr}^{2}}-\dfrac{{\mathcal{S}}(v)}{\hat{u}_{th}^{2}}$
(64o) $\displaystyle\Gamma_{-1b}(v)$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle-\dfrac{A_{r}^{3}}{U_{b}^{2}}{\mathcal{S}}\left(v+\dfrac{Q}{A_{r}}\right)-\dfrac{{\mathcal{S}}(v)}{\hat{u}_{pr}^{2}}+F_{hM}$
(64p) $\displaystyle\Gamma_{-2a}(v)$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$ $\displaystyle-
A_{r}P_{M}-\dfrac{{\mathcal{S}}(v)}{\hat{u}_{pr}^{2}}-\dfrac{{\mathcal{S}}(v)}{\hat{u}_{th}^{2}}$
(64q) $\displaystyle\Gamma_{-2b}(v)$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$ $\displaystyle-
A_{r}P_{M}-\dfrac{{\mathcal{S}}(v)}{\hat{u}_{pr}^{2}}+F_{hM}$ (64r)
$\displaystyle\Gamma_{-3}(v)$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle-\dfrac{{\mathcal{S}}(v)}{\hat{u}_{th}^{2}}.$ (64s)

Here, the divisions by $\hat{u}_{\ast}^{2}$ or $U_{b}^{2}$ do not cause
troubles even when they are zeros because (64c) implicitly includes upper and
lower bounds. One workaround in the implementation may be to replace the
divisions by $\hat{u}_{\ast}^{2}$ by those by
$\max(\varepsilon,\hat{u}_{\ast}^{2})$ where $\varepsilon$ is a small positive
number close to the machine epsilon.

Table 1: States of valves at each segment of $\Gamma(v)$. hR = head-side relief valve; hSC = head-side suction check valve; pR = pump relief valve; rSC = rod-side suction relief valve; rR = rod-side relief valve; o = open; x = closed; - = either open or closed. $\Gamma$ | $v$ | $\hat{{\boldsymbol{u}}}$ | hR | hSC | pR | rSC | rR
---|---|---|---|---|---|---|---
$-F_{rM}$ | $+$ | ${\boldsymbol{\mathcal{U}}}_{+}\cup{\boldsymbol{\mathcal{U}}}_{0}\cup{\boldsymbol{\mathcal{U}}}_{-}$ | x | o | - | x | o
$\Gamma_{+3}$ | $+$ | ${\boldsymbol{\mathcal{U}}}_{+}$ | x | o | - | x | x
$\Gamma_{+2b}$ | $+$ | ${\boldsymbol{\mathcal{U}}}_{+}$ | x | x | o | x | o
$\Gamma_{+2a}$ | $+$ | ${\boldsymbol{\mathcal{U}}}_{+}$ | x | x | o | x | x
$\Gamma_{+1b}$ | $+$ | ${\boldsymbol{\mathcal{U}}}_{+}$ | x | x | x | x | o
$\Gamma_{+1a}$ | $+$ | ${\boldsymbol{\mathcal{U}}}_{+}$ | x | x | x | x | x
$\Gamma_{+0b}$ | $+$ | ${\boldsymbol{\mathcal{U}}}_{+}$ | o | x | - | x | o
$\Gamma_{+0a}$ | $+$ | ${\boldsymbol{\mathcal{U}}}_{+}$ | o | x | - | x | x
$[\Gamma_{+}(0),F_{hM}]$ | $0$ | ${\boldsymbol{\mathcal{U}}}_{+}$ | x | x | - | x | x
$[-F_{rM},F_{hM}]$ | $0$ | ${\boldsymbol{\mathcal{U}}}_{0}$ | x | x | - | x | x
$[-F_{rM},\Gamma_{-}(0)]$ | $0$ | ${\boldsymbol{\mathcal{U}}}_{-}$ | x | x | - | x | x
$\Gamma_{-0a}$ | $-$ | ${\boldsymbol{\mathcal{U}}}_{-}$ | x | x | - | x | o
$\Gamma_{-0b}$ | $-$ | ${\boldsymbol{\mathcal{U}}}_{-}$ | o | x | - | x | o
$\Gamma_{-1a}$ | $-$ | ${\boldsymbol{\mathcal{U}}}_{-}$ | x | x | x | x | x
$\Gamma_{-1b}$ | $-$ | ${\boldsymbol{\mathcal{U}}}_{-}$ | o | x | x | x | x
$\Gamma_{-2a}$ | $-$ | ${\boldsymbol{\mathcal{U}}}_{-}$ | x | x | o | x | x
$\Gamma_{-2b}$ | $-$ | ${\boldsymbol{\mathcal{U}}}_{-}$ | o | x | o | x | x
$\Gamma_{-3}$ | $-$ | ${\boldsymbol{\mathcal{U}}}_{-}$ | x | x | - | o | x
$F_{hM}$ | $-$ | ${\boldsymbol{\mathcal{U}}}_{+}\cup{\boldsymbol{\mathcal{U}}}_{0}\cup{\boldsymbol{\mathcal{U}}}_{-}$ | o | x | - | o | x

It should be noted that, at $v=0$, the function $\Gamma(v)$ is set-valued and
its value is the closed set $[\Gamma_{+}(0),\Gamma_{-}(0)]$. The boundaries
$\Gamma_{+}(0)$ and $\Gamma_{-}(0)$ are obtained from a straightforward
derivation as

$\displaystyle\Gamma_{+}(0)=\Gamma_{h+}(0)-\Gamma_{r+}(0),\quad\Gamma_{-}(0)=\Gamma_{h-}(0)-\Gamma_{r-}(0)$
(65a) where $\displaystyle\Gamma_{h+}(0)$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}\min(F_{hM},A_{h}P_{M},A_{h}Q^{2}/U_{b}^{2})&\mbox{if
}\hat{u}_{ph}>0\\\ 0&\mbox{if }\hat{u}_{ph}=0\\\ \end{array}\right.$ (65d)
$\displaystyle\Gamma_{r+}(0)$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}0&\mbox{if }\hat{u}_{tr}>0\\\
F_{rM}&\mbox{if }\hat{u}_{tr}=0\\\ \end{array}\right.$ (65g)
$\displaystyle\Gamma_{h-}(0)$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}0&\mbox{if }\hat{u}_{th}>0\\\
F_{hM}&\mbox{if }\hat{u}_{th}=0\\\ \end{array}\right.$ (65j)
$\displaystyle\Gamma_{r-}(0)$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}\min(F_{rM},A_{r}P_{M},A_{r}Q^{2}/U_{b}^{2})&\mbox{if
}\hat{u}_{pr}>0\\\ 0&\mbox{if }\hat{u}_{pr}=0.\\\ \end{array}\right.$ (65m)

As can be seen from the expression (64), the function $\Gamma(v)$ is composed
of many segments. Interestingly, each segment has a clear physical
interpretation. For example, in the segments listed in (64f)-(64l), the last
term, either $-F_{rM}$ or $-{\mathcal{S}}(v)/\hat{u}_{tr}^{2}$, represents the
state of the rod-side relief valve, either closed or open, respectively. In
this way of consideration, we can summarize the valve states at each curve
segment as in Table 1.

###  3.4 Numerical Examples

Figure 3: Numerical examples of the quasistatic map $f\in\Gamma(v)$. (a)
$u_{c}=0.5$, $u_{b}=0.2$; (b) $u_{c}=-0.5$, $u_{b}=0.2$; (c) $u_{c}\in[-1,1]$,
$u_{b}=0.2$; (d) $u_{c}=0.7$, $u_{b}\in[0,1]$. The velocity $v$ is positive
when the rod is extending. The force $f$ is positive when the external force
is compressive and the actuator force acts to extend the rod.

Some numerical examples are now presented. We consider an asymmetric hydraulic
cylinder with the following parameter values:

$\displaystyle C_{*}=0.6,\ a_{*}=0.0001~{}\mbox{m${}^{2}$},\
\rho=850~{}\mbox{kg/m${}^{3}$},$ $\displaystyle
A_{r}=0.012~{}\mbox{m${}^{2}$},\ A_{h}=0.024~{}\mbox{m${}^{2}$},$
$\displaystyle P_{rM}=40~{}\mbox{MPa},\ P_{hM}=42~{}\mbox{MPa},\
P_{M}=36~{}\mbox{MPa},$ $\displaystyle
Q=500~{}\mbox{L/min}=0.00833~{}\mbox{m${}^{3}$/s}.$ (66)

A common control command $u_{c}\in[-1,1]$ is used to set the openings of the
main control valves as follows:

$\displaystyle u_{ph}=u_{tr}=\max(u_{c},0),\ u_{pr}=u_{th}=\max(-u_{c},0).$
(67)

The opening of the bleed valve was fixed at $u_{b}=0.2$ unless otherwise
noted.

Results with different values of $u_{c}$ and $u_{b}$ are shown in Fig. 3. Fig.
3(a) shows the function $\Gamma$ and its segments with a positive $u_{c}$,
while Fig. 3(b) shows those with a negative $u_{c}$. They show that
$\Gamma(v)$ is always a decreasing function of $v$ and it is set-valued at
$v=0$. Fig. 3(c) shows how the function $\Gamma(v)$ varies according to the
change in $u_{c}$. It shows that, at a constant external force, the velocity
$v$ increases as $u_{c}$ increases, which is consistent with the behavior of
real hydraulic actuator. Fig. 3(d) shows that the velocity $v$ decreases as
the bleed valve opening $u_{b}$ increases. It can be explained by the fact
that the oil flow into the actuator decreases as the bleed valve opens.

##  4 Incorporation into Multibody Simulators

###  4.1 Virtual Viscoelastic Element

The previous section has modeled a hydraulic actuator as a set-valued function
$\Gamma$ from the velocity $v$ to the force $f$. Because of the set-
valuedness, the function $\Gamma$ is not convenient for the use in
simulations. The set-valuedness at $v=0$ is an important feature that cannot
be neglected because the actuator actually stops when the oil flow is blocked
by the valves. This section shows an approach to incorporate the set-valued
function $\Gamma$ in multibody simulators.

Figure 4: Differential-algebraic relaxation using a virtual viscoelastic
element. (a) The original system with the algebraic constraint
$f\in\Gamma(\dot{\ell})$. (b) The approximate (relaxed) system with the
differential-algebraic constraint (68).

The quasistatic model presented in the previous section describes the system
of Fig. 4(a) as an algebraic constraint $f\in\Gamma(\dot{p})$ where $p$ is the
displacement of the actuator rod and $v=\dot{p}$. To deal with the set-
valuedness, we consider installing a virtual viscoelastic element connected
with the actuator as illustrated in Fig. 4(b). Let $\ell$ be the total
displacement including the viscoelastic element and $p$ be the displacement of
the actuator rod. The stretch of the viscoelastic element is $\ell-p$. One can
see that, if the viscoelastic element is stiff enough, $p\approx\ell$ is
maintained and thus the system of Fig. 4(b) approximates the system of Fig.
4(a). The dynamics of the system of Fig. 4(b) can be described as follows:

$\displaystyle f$ $\displaystyle=$ $\displaystyle
K(p-\ell)+B(\dot{p}-\dot{\ell})$ (68a) $\displaystyle f$ $\displaystyle\in$
$\displaystyle\Gamma(\dot{p}).$ (68b)

Here, $K$ and $B$ are the stiffness and the viscosity of the virtual
viscoelastic element. Because there is no inertia between the actuator and the
viscoelastic element, $f$ in (68a), determined by the viscoelastic element, is
always the same as $f$ in (68b), acting on the actuator. Although some
discussions on the choice of $K$ and $B$ will be given in Remark 2, a simple
guideline for a better approximation (i.e., $p\approx\ell$) is that the
stiffness $K$ should be as high as the simulation is stable, and that the
viscosity $B$ should be high enough to damp the oscillation in the relative
displacement $p-\ell$. The expression (68) can be seen as a differential-
algebraic inclusion (DAI) with respect to $p$. The quantities $\ell$ and
$\dot{\ell}$ are given as inputs, and $\dot{p}$ and $f$ are determined by the
pair (68) of algebraic constraints.

To provide the solution of these algebraic constraints (68), let us consider a
function $\Lambda$ that satisfies the following relation:

$\displaystyle\beta v+\bar{f}\in\Gamma(v)$ $\displaystyle\iff$ $\displaystyle
v=\Lambda(\beta,\bar{f})$ (69)

where $\beta>0$. Because $\Gamma(v)$ is a decreasing function of $v$ and
$\beta v+\bar{f}$ is a strictly increasing function of $v$, the algebraic
inclusion $\beta v+\bar{f}\in\Gamma(v)$ always has a unique solution with
respect to $v$, and thus $\Lambda(\beta,\bar{f})$ is a globally single-valued
function. Fortunately, $\Lambda(\beta,\bar{f})$ can be obtained analytically
as will be shown in the next Section 4.2. Considering that (68a) is equivalent
to $f=B\dot{p}+\bar{f}$ where $\bar{f}=K(p-\ell)-B\dot{\ell}$, the DAI (68)
can be equivalently rewritten as follows:

$\displaystyle\dot{p}$ $\displaystyle=$
$\displaystyle\Lambda(B,K(p-\ell)-B\dot{\ell})$ (70a) $\displaystyle f$
$\displaystyle=$ $\displaystyle K(p-\ell)+B(\dot{p}-\dot{\ell}).$ (70b)

The expression (70) is algebraically equivalent to the expression (68), and
can be seen as a closed-form solution of the algebraic problem (68). Moreover,
because $\Lambda$ is a continuous, single-valued function, (70) is only an
ordinary differential equation (ODE) representing a first-order system of
which the state is $p$. That is, (70) can be used for simulations combined
with common numerical integration schemes.

Another numerical scheme for the DAI (68) can be obtained through the
discretization along time. With the implicit (backward) Euler discretization,
i.e., $v=\dot{p}\approx(p_{k}-p_{k-1})/h$ and
$\dot{\ell}\approx(\ell_{k}-\ell_{k-1})/h$, one can discretize (68) as
follows:

$\displaystyle f_{k}$ $\displaystyle=$
$\displaystyle(B+hK)v+K(p_{k-1}-\ell_{k})-B(\ell_{k}-\ell_{k-1})/h$ (71a)
$\displaystyle f_{k}$ $\displaystyle\in$ $\displaystyle\Gamma(v)$ (71b)
$\displaystyle p_{k}$ $\displaystyle=$ $\displaystyle p_{k-1}+hv$ (71c)

where $k$ denote the discrete-time index and $h$ denotes the timestep size.
The pair of constraints (71a) and (71b) poses an algebraic problem with
respect to $v$ and $f_{k}$, and the function $\Lambda$ provides the solutions.
That is, we can solve (71) in the following algorithm:

$\displaystyle\bar{f}$ $\displaystyle:=$ $\displaystyle
K(p_{k-1}-\ell_{k})-B(\ell_{k}-\ell_{k-1})/h$ (72a) $\displaystyle v$
$\displaystyle:=$ $\displaystyle\Lambda(B+hK,\bar{f})$ (72b) $\displaystyle
p_{k}$ $\displaystyle:=$ $\displaystyle p_{k-1}+hv$ (72c) $\displaystyle
f_{k}$ $\displaystyle:=$ $\displaystyle\bar{f}+(B+hK)v.$ (72d)

This algorithm accepts the inputs $\\{\ell_{k},p_{k-1},\ell_{k-1}\\}$ and
provides the output $f_{k}$, and can be seen as a numerical integration scheme
for the DAI (68). In the simulation, the values of $\\{p_{k},\ell_{k}\\}$
should be stored and used as $\\{p_{k-1},\ell_{k-1}\\}$ in the next timestep.

In summary, the original set-valued constraint $f\in\Gamma(\dot{\ell})$ is
approximated by the DAI (68), which is implementable as the ODE (70) or the
discrete-time algorithm (72). For the implementation, the continuous-time form
(70) would be preferred for rather complicated simulations with sophisticated
ODE solvers, while the Euler-based discrete-time form (72) would be convenient
for rather simple simulations that do not require high accuracy. This
approximation scheme, which relaxes a set-valued algebraic constraint by a
DAI, is what Kikuuwe [15] has called a differential algebraic relaxation. The
scheme has been utilized to deal with Coulomb friction in simulators [17, 31,
30, 16] and to implement sliding mode controllers to discrete-time systems
[18, 14].

######  Remark 2.

The choice of $K$ and $B$ can be discussed from two different points of view.
On the one hand, if one needs to approximate the ideal situation where the
actuator is rigidly connected to external components (such as links of an
excavator), the stiffness $K$ should be set as high as the simulation is
stable. On the other hand, if one needs to simulate a particular real system
that comprises the compliance of the components and the compressibility of the
oil, the values of $K$ and $B$ should be chosen so that the simulator’s
response is close to that of the real system. In this case, the values of $K$
and $B$ may be considered as the aggregation of the compliance and damping of
the components and the oil. The tuning procedure might be performed by
explicitly considering these factors, or on a trial-and-error basis, comparing
the data from the simulator and the real system.

###  4.2 Analytical Form of Function $\Lambda$

Now we present an analytical form of $\Lambda(\beta,\bar{f})$, of which the
definition has been given only implicitly as in (69). Let $v_{\ast}$ denote
the solution of $\beta v+\bar{f}=\Gamma_{\ast}(v)$ with respect to $v$ where
$\Gamma_{\ast}(v)$ are those listed in (64). By observing the definitions of
$\Gamma_{\ast}(v)$ and using Proposition 1 and Proposition 2, we can obtain
the followings:

$\displaystyle v_{+0a}$ $\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle\Phi_{A}\left(\beta,\bar{f}-F_{hM},\hat{u}_{tr}^{2}\right)$
(73a) $\displaystyle v_{+0b}$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle(F_{hM}-F_{rM}-\bar{f})/\beta$ (73b) $\displaystyle v_{+1a}$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle\Phi_{B}\left(\beta,\bar{f},\psi(\hat{u}_{ph},\hat{u}_{tr}),\dfrac{U_{b}^{2}}{A_{h}^{3}},\dfrac{Q}{A_{h}}\right)$
(73c) $\displaystyle v_{+1b}$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle\Phi_{B}\left(\beta,\bar{f}+F_{rM},\hat{u}_{ph}^{2},\dfrac{U_{b}^{2}}{A_{h}^{3}},\dfrac{Q}{A_{h}}\right)$
(73d) $\displaystyle v_{+2a}$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle\Phi_{A}\left(\beta,\bar{f}-A_{h}P_{M},\psi(\hat{u}_{ph},\hat{u}_{tr})\right)$
(73e) $\displaystyle v_{+2b}$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle\Phi_{A}\left(\beta,\bar{f}-A_{h}P_{M}+F_{rM},\hat{u}_{ph}^{2}\right)$
(73f) $\displaystyle v_{+3}$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle\Phi_{A}\left(\beta,\bar{f},\hat{u}_{tr}^{2}\right)$ (73g)
$\displaystyle v_{rM}$ $\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle(-F_{rM}-\bar{f})/\beta$ (73h) $\displaystyle v_{-0a}$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle\Phi_{A}\left(\beta,\bar{f}+F_{rM},\hat{u}_{th}^{2}\right)$
(73i) $\displaystyle v_{-0b}$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle(F_{hM}-F_{rM}-\bar{f})/\beta$ (73j) $\displaystyle v_{-1a}$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle\Phi_{B}\left(\beta,\bar{f},\psi(\hat{u}_{pr},\hat{u}_{th}),\dfrac{U_{b}^{2}}{A_{r}^{3}},-\dfrac{Q}{A_{r}}\right)$
(73k) $\displaystyle v_{-1b}$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle\Phi_{B}\left(\beta,\bar{f}-F_{hM},\hat{u}_{pr}^{2},\dfrac{U_{b}^{2}}{A_{r}^{3}},-\dfrac{Q}{A_{r}}\right)$
(73l) $\displaystyle v_{-2a}$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle\Phi_{A}\left(\beta,\bar{f}+A_{r}P_{M},\psi(\hat{u}_{pr},\hat{u}_{th})\right)$
(73m) $\displaystyle v_{-2b}$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle\Phi_{A}\left(\beta,\bar{f}+A_{r}P_{M}-F_{hM},\hat{u}_{pr}^{2}\right)$
(73n) $\displaystyle v_{-3}$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle\Phi_{A}\left(\beta,\bar{f},\hat{u}_{th}^{2}\right)$ (73o)
$\displaystyle v_{hM}$ $\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle(F_{hM}-\bar{f})/\beta.$ (73p)

Then, by using Theorem 2 and Theorem 3, the function $\Lambda$ can be obtained
as follows:

$\displaystyle\Lambda(\beta,\bar{f})$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle\left\\{\begin{array}[]{ll}v_{hM}&\mbox{if
}\bar{f}>F_{hM}\,\wedge\,\hat{{\boldsymbol{u}}}\in({\boldsymbol{\mathcal{U}}}_{+}\cup\,{\boldsymbol{\mathcal{U}}}_{0})\\\
\Lambda_{+}(\beta,\bar{f})&\mbox{if
}\bar{f}<\Gamma_{+}(0)\,\wedge\,\hat{{\boldsymbol{u}}}\in{\boldsymbol{\mathcal{U}}}_{+}\\\
v_{rM}&\mbox{if
}\bar{f}<-F_{rM}\,\wedge\,\hat{{\boldsymbol{u}}}\in({\boldsymbol{\mathcal{U}}}_{-}\cup\,{\boldsymbol{\mathcal{U}}}_{0})\\\
\Lambda_{-}(\beta,\bar{f})&\mbox{if
}\bar{f}>\Gamma_{-}(0)\,\wedge\,\hat{{\boldsymbol{u}}}\in{\boldsymbol{\mathcal{U}}}_{-}\\\
0&\mbox{if }\Gamma_{+}(0)\leq\bar{f}\leq\Gamma_{-}(0)\end{array}\right.$ (74f)
where $\displaystyle\Lambda_{+}(\beta,\bar{f})$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle\max(\min(\max(v_{+0a},v_{+0b}),\max(v_{+1a},v_{+1b}),\max(v_{+2a},v_{+2b})),v_{+3},v_{rM})$
(74g) $\displaystyle\Lambda_{-}(\beta,\bar{f})$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle\min(\max(\min(v_{-0a},v_{-0b}),\min(v_{-1a},v_{-1b}),\min(v_{-2a},v_{-2b})),v_{-3},v_{hM}).$
(74h)

Here, $\Gamma_{+}(0)$ and $\Gamma_{-}(0)$ are those given by (65).

###  4.3 Numerical Examples

Figure 5: A single-DOF arm driven by a hydraulic actuator.

Figure 6: Inputs (a) and results (b)-(f) of the simulation of the system of
Fig. 6. The actuator force $f$ is positive when it acts to extend the rod. The
actuator force $f$ saturates at $t=7.404$ s, indicated by the vertical lines.

The method is illustrated with a simple dynamics simulation of a single
degree-of-freedom (DOF) arm system shown in Fig. 6. The equation of motion is
written as follows:

$\displaystyle(J+ML_{g}^{2})\ddot{\theta}$ $\displaystyle=$
$\displaystyle{\boldsymbol{r}}_{g}\times\left[\begin{array}[]{c}0\\\
-Mg\end{array}\right]+{\boldsymbol{r}}_{m}\times\left(\dfrac{{\boldsymbol{r}}_{m}-{\boldsymbol{r}}_{b}}{\|{\boldsymbol{r}}_{m}-{\boldsymbol{r}}_{b}\|}f\right)+{\boldsymbol{r}}_{f}\times\left[\begin{array}[]{c}0\\\
f_{ey}\end{array}\right]$ (79) $\displaystyle\ell$ $\displaystyle=$
$\displaystyle\|{\boldsymbol{r}}_{m}-{\boldsymbol{r}}_{b}\|$ (80)
$\displaystyle f$ $\displaystyle=$ $\displaystyle
K(\ell-p)+B(\dot{\ell}-\dot{p})$ (81) $\displaystyle f$ $\displaystyle\in$
$\displaystyle\Gamma(\dot{p})$ (82)

where ${\boldsymbol{r}}_{*}$ are those indicated in Fig. 6, which are
calculaed as follows:

$\displaystyle{\boldsymbol{r}}_{g}=\left[\begin{array}[]{c}L_{g}\cos\theta\\\
L_{g}\sin\theta\end{array}\right],\quad{\boldsymbol{r}}_{m}=\left[\begin{array}[]{c}-L_{m}\cos(\theta-\alpha)\\\
-L_{m}\sin(\theta-\alpha)\end{array}\right],\quad{\boldsymbol{r}}_{f}=\left[\begin{array}[]{c}L_{f}\cos\theta\\\
L_{f}\sin\theta\end{array}\right].$ (89)

Here, $\times$ stands for a pseudo cross-product of two-dimensional vectors,
$[a_{1},a_{2}]^{T}\times[b_{1},b_{2}]^{T}=a_{1}b_{2}-a_{2}b_{1}$. The
parameters were set as $L_{g}=1.5$ m, $L_{m}=0.6$ m, $\alpha=\pi/4$, $L_{f}=3$
m, $M=2000$ kg, and $J=5000$ kg$\cdot$m2. The parameters for the virtual
viscoelastic elements were chosen as $K=5\times 10^{7}$ N/m and $B=2.5\times
10^{6}$ N/m, which as high as the simulation was numerically stable. The
parameters of the actuator were set the same as in Section 3.4. The openings
of the main control valves were given by (67) with the command
$u_{c}\in[-1,1]$, and the bleed valve opening was set as $u_{b}=0.3$. The
timestep size was set as $h=0.001$ s. At every timestep, the following
procedure was performed:

  1. 1.

Calculate ${\boldsymbol{r}}_{g}$, ${\boldsymbol{r}}_{m}$, and
${\boldsymbol{r}}_{f}$ according to (89) with $\theta$ determined by the
previous timestep, and calculate $\ell$ according to (80).

  2. 2.

Calculate $f$ and update $p$ according to the algorithm (72) with the command
$u_{c}$ and the obtained $\ell$.

  3. 3.

Update $\theta$ and $\dot{\theta}$ according to the equation of motion (79)
with the inputs $f$ and $f_{ey}$.

The initial value of $p$ was set equal to that of $\ell$.

In the simulation, the command $u_{c}$ and $f_{ey}$ were given as indicated in
Fig. 6(a). The results are shown in Figs. 6(b)-(f). From $t=0$ s to $3$ s, the
valves are closed and thus the arm holds its angle even under the fluctuating
external force $f_{ey}$. From $t=3$ s to $6$ s, the arm moves up due the
negative $u_{c}$ (to retract the rod). Fig. 6(c) shows that the velocity is
more fluctuating in this period than when $u_{c}=0$. This behavior is
consistent with the fact that hydraulic actuators become more compliant to
external forces when the valves are open and the oil is permitted to flow.
After $t=6$ s, the valves are closed again, and the increasing external force
does not cause movement of the arm until $t=7.404$ s, with the increasing
actuator force $f$ balancing the external force. At $t=7.404$ s, the actuator
force $f$ reaches the value of $-F_{rM}$ at which the rod-side relief valve
opens, and the arm starts to move down. In the whole period, the displacement
$p-\ell$ of the viscoelastic element is maintained small (less than a few
centimeters), exhibiting the validity of the relaxation scheme introduced in
Section 4.1. In conclusion, the results support the usefulness of the proposed
model and the relaxation scheme for reproducing qualitative features of
hydraulic actuators.

######  Remark 3.

The transient responses seen in Fig. 6, such as the spikes in the actuator
force $f$ at $t\approx 3$ s and $t\approx 6$ s, are determined by the values
of $K$ and $B$. As discussed in Remark 2, one needs to tune the values of $K$
and $B$ so that the simulator exhibits similar responses to the real system.

##  5 Extension 1: Regeneration Circuit

###  5.1 Quasistatic Model

Figure 7: Regeneration pipeline added to the circuit of Fig. 1.

In some practical applications, the circuit of Fig. 1 includes an additional
pathway from the rod-side chamber to the head-side chamber to make the
extending movement faster. This section considers the circuit of Fig. 7, which
is an extension of the circuit of Fig. 1. We assume that the regeneration
pipeline has a valve that flows only from the head side to the rod side, and
assume that the cross-sectional areas of the cylinder satisfy

$\displaystyle A_{h}\geq A_{r}.$ (90)

There must be some cylinders that do not satisfy (90), but the removal of this
assumption is left for future study. Let $a_{a}$ is the maximum opening area
(m2) of the regeneration valve, and let $u_{a}\in[0,1]$ be the dimensionless
input value, which is the ratio of the valve opening area to its maximum value
$a_{a}$. Then, the flowrate of the oil through the regeneration pipeline is
written as follows:

$\displaystyle Q_{a}$ $\displaystyle=$ $\displaystyle
c_{a}u_{a}\max({\mathcal{R}}(F_{r}/A_{r}-F_{h}/A_{h},0))$ (91)

where $c_{a}\stackrel{{\scriptstyle\Delta}}{{=}}C_{a}a_{a}\sqrt{2/\rho}$ where
$C_{a}$ is the discharge coefficient of the regeneration valve, typically
around $0.6$ or $0.7$.

With this regeneration pipeline and its flowrate $Q_{a}$, (47) is extended
into the following:

$\displaystyle-v+\hat{u}_{ph}{\mathcal{R}}(A_{h}P_{c}-F_{h})-\hat{u}_{th}{\mathcal{R}}(F_{h})\in{\mathcal{N}}_{[0,F_{hM}]}(F_{h})-Q_{a}/A_{h}$
(92a) $\displaystyle
v+\hat{u}_{pr}{\mathcal{R}}(A_{r}P_{c}-F_{r})-\hat{u}_{tr}{\mathcal{R}}(F_{r})\in{\mathcal{N}}_{[0,F_{rM}]}(F_{r})+Q_{a}/A_{r}$
(92b) $\displaystyle
A_{h}\hat{u}_{ph}{\mathcal{R}}(A_{h}P_{c}-F_{h})+A_{r}\hat{u}_{pr}{\mathcal{R}}(A_{r}P_{c}-F_{r})\in{\mathcal{N}}_{(-\infty,P_{c}]}(P)$
(92c) $\displaystyle Q\in
U_{b}{\mathcal{R}}(P)+A_{h}\hat{u}_{ph}{\mathcal{R}}(A_{h}P_{c}-F_{h})+A_{r}\hat{u}_{pr}{\mathcal{R}}(A_{r}P_{c}-F_{r})+{\mathcal{N}}_{(-\infty,P_{M}]}(P)$
(92d) $\displaystyle f=F_{h}-F_{r}.$ (92e)

With $Q_{a}=0$, which results from $u_{a}=0$, (92) reduces to (47). Let us
assume that $u_{a}$ is set positive only when
$\hat{{\boldsymbol{u}}}\in{\boldsymbol{\mathcal{U}}}_{+}$ because it is the
case in many practical hydraulic circuits. In addition, we assume that
$u_{tr}>0$ for the simplicity. Under these conditions, (92) and (91) reduce to
the following:

$\displaystyle\hat{A}v_{a}-v+\hat{u}_{ph}{\mathcal{R}}(A_{h}P_{c}-F_{h})\in{\mathcal{N}}_{[0,F_{hM}]}(F_{h})$
(93a) $\displaystyle
v-v_{a}-\hat{u}_{tr}{\mathcal{R}}(F_{r})\in{\mathcal{N}}_{[0,F_{rM}]}(F_{r})$
(93b) $\displaystyle
A_{h}\hat{u}_{ph}{\mathcal{R}}(A_{h}P_{c}-F_{h})\in{\mathcal{N}}_{(-\infty,P_{c}]}(P)$
(93c) $\displaystyle Q\in
U_{b}{\mathcal{R}}(P)+A_{h}\hat{u}_{ph}{\mathcal{R}}(A_{h}P_{c}-F_{h})+{\mathcal{N}}_{(-\infty,P_{M}]}(P)$
(93d) $\displaystyle f=F_{h}-F_{r}$ (93e) $\displaystyle
v_{a}=\hat{u}_{a}\max({\mathcal{R}}(F_{r}-\hat{A}F_{h},0))$ (93f)

where $v_{a}\stackrel{{\scriptstyle\Delta}}{{=}}Q_{a}/A_{r}$,
$\hat{u}_{a}\stackrel{{\scriptstyle\Delta}}{{=}}c_{a}u_{a}/A_{r}^{3/2}$ and
$\hat{A}=A_{r}/A_{h}$.

By carefully observing (93), one can see that $v_{a}>0$ implies $F_{r}>0$ from
(93f), which implies $v_{a}\leq v$ from (93b). Its contraposition is that
$v_{a}>v$ implies $v_{a}=0$. Therefore, (93) imposes the condition $0\leq
v_{a}\leq\max(0,v)$. When $v<0$, the solutions are obviously $v_{a}=0$ and
$f\in\Gamma(v)$. Therefore, hereafter we consider only the case $v>0$, in
which $v_{a}>0$ may take place. By using the functions defined in (61) and
(62), the first four equations of (93) can be rewritten as
$F_{h}\in\Gamma_{h}(v-\hat{A}v_{a})$ and $F_{r}\in\Gamma_{r}(v-v_{a})$ with
$P_{c}$ and $P$ being eliminated. Considering the definitions of $\Gamma_{h}$
and $\Gamma_{r}$, because of the conditions $v-v_{a}\geq 0$, $v>0$, $u_{tr}>0$
and $\hat{A}\in(0,1)$, $\Gamma_{h}$ and $\Gamma_{r}$ are always single-valued
and thus can be replaced by $\Gamma_{h+}$ and $\Gamma_{r+}$, respectively.
Therefore, (93) can be rewritten as follows:

$\displaystyle\Xi_{v}(v,v_{a})\in-{\mathcal{N}}_{[0,\infty)}(v_{a})$ (94a)
$\displaystyle f=\Gamma_{h+}(v-\hat{A}v_{a})-\Gamma_{r}(v-v_{a})$ (94b)

where

$\displaystyle\Xi_{v}(v,v_{a})$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle{\mathcal{S}}(v_{a})-\hat{u}_{a}^{2}(\Gamma_{r+}(v-v_{a})-\hat{A}\Gamma_{h+}(v-\hat{A}v_{a})).$
(95)

This expression can be seen as an algebraic problem regarding $\\{f,v_{a}\\}$
with a given $v$.

The function $\Xi_{v}(v,v_{a})$ is an increasing function of $v_{a}$ and it
satisfies $\Xi_{v}(v,v)>0$ because of $\Gamma_{r+}(0)=0$. With the algebraic
constraint (94a), $\Xi_{v}(v,0)\geq 0$ implies that the solution is $v_{a}=0$,
Meanwhile, if $\Xi_{v}(v,0)<0$, the solution can be found within the region
$v_{a}\in[0,v]$ by simple root-finding schemes. Once the solution $v_{a}$ is
obtained, $f$ is obtained by (94b).

In conclusions, $f$ and $v_{a}$ satisfying (94) are obtained by the following
functions:

$\displaystyle f\in\Gamma_{\mathrm{reg}}(v)$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle\left\\{\begin{array}[]{lr}\Gamma(v)&\mbox{if
}u_{a}=0\,\vee\,\hat{{\boldsymbol{u}}}\not\in{\boldsymbol{\mathcal{U}}}_{+}\,\vee\,v\leq
0\,\vee\,\Xi_{v}(v,0)\geq 0\\\
\Gamma_{h+}(v-\hat{A}\widehat{v}_{a}(v))-\Gamma_{r+}(v-\widehat{v}_{a}(v))&\mbox{otherwise}\\\
\end{array}\right.$ (98) $\displaystyle v_{a}=\widehat{v}_{a}(v)$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle\left\\{\begin{array}[]{lr}0&\mbox{if
}u_{a}=0\,\vee\,\hat{{\boldsymbol{u}}}\not\in{\boldsymbol{\mathcal{U}}}_{+}\,\vee\,v\leq
0\,\vee\,\Xi_{v}(v,0)\geq 0\\\
\mbox{FindRoot}(\Xi_{v}(v,\bullet),[0,v])&\mbox{otherwise.}\end{array}\right.$
(101)

Here, “FindRoot” is a function that finds a root of the argument function
within the range specified by the second argument. This computation can be
performed with a common iterative method, such as the Bisection method or the
false position method, because $\Xi_{v}(v,v_{a})$ is continuous and monotonic.
It may also be possible to find analytical methods because $\Xi_{v}(v,v_{a})$
is only a piece-wise parabolic function. The function $\Gamma_{\mathrm{reg}}$
can be seen as an extension of the quasistatic map $\Gamma$ given in Section
3.3.

Figure 8: Numerical examples of the quasistatic map
$f\in\Gamma_{\mathrm{reg}}(v)$ and $v_{a}=\widehat{v}_{a}(v)$ of the circuit
of Fig. 7, which includes a regeneration pipeline. The parameters are set the
same as in Section 3.4 except $C_{a}=0.6$ and $a_{a}=0.0001$ m2. The valve
openings are set as $u_{c}=0.5$ (see (67)) and $u_{b}=0.2$.

Some numerical examples of $\Gamma_{\mathrm{reg}}$ and $\widehat{v}_{a}$ are
presented in Fig. 8. It can be seen that, as $u_{a}$ increases (i.e., as the
regeneration valve opens), the extending velocity of the rod tends to increase
especially in the region $f<0$, i.e., under stretching external forces. It can
also be seen that the flowrate $Q_{a}$ through the regeneration valve
increases as $u_{a}$ and $v$ increase. These features are consistent with
those of actual hydraulic actuators.

###  5.2 Incorporation into Simulators

As has been discussed in Section 4, multibody simulations involving the
quasistatic map $f\in\Gamma(v)$ require its correspondent $\Lambda$ defined by
(69). In the same manner, now we consider obtaining the function
$\Lambda_{\mathrm{reg}}$ that satisfies the following:

$\displaystyle\beta v+\bar{f}\in\Gamma_{\mathrm{reg}}(v)$ $\displaystyle\iff$
$\displaystyle v=\Lambda_{\mathrm{reg}}(\beta,\bar{f}).$ (102)

To this end, we consider the algebraic problem (94) with $f$ being replaced by
$f=\beta v+\bar{f}$ and $v$ being treated as an unknown, which is written as
follows:

$\displaystyle\Xi_{v}(v,v_{a})$ $\displaystyle\in$
$\displaystyle-{\mathcal{N}}_{[0,\infty)}(v_{a})$ (103a)
$\displaystyle\Xi_{f}(v,v_{a})$ $\displaystyle=$ $\displaystyle 0$ (103b)

where

$\displaystyle\Xi_{f}(v,v_{a})$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$ $\displaystyle\beta
v+\bar{f}-\Gamma_{h+}(v-\hat{A}v_{a})+\Gamma_{r+}(v-v_{a}).$ (104)

The expression (103) is an algebraic problem regarding $\\{v,v_{a}\\}$ with a
given $\bar{f}$.

The problem (103) is illustrated in Fig. 10, in which the sets of
$\\{v,v_{a}\\}$ satisfying (103a) and (103b) are shown as curves. The
functions $\Xi_{v}(v_{a},v)$ and $\Xi_{f}(v_{a},v)$ are increasing with
respect to $v_{a}$ and $v$, respectively. The region $0\leq v_{a}<v$ should be
searched for the solution $\\{v,v_{a}\\}$. Because $f=\beta
v+\bar{f}\in[-F_{rM},F_{hM}]$, one can see that the solution exists in the
region $(-F_{rM}-\bar{f})/\beta<v<(F_{hM}-\bar{f})/\beta$. In addition, the
solution $v$ of $\Xi_{f}(v,v_{a})=0$ with $v_{a}>0$ is larger than
$\Lambda(\beta,\bar{f})$, which is the solution with $v_{a}=0$. Therefore, the
trapezoidal area in Fig. 10 is the region that must be searched for the
solution $\\{v,v_{a}\\}$.

As shown in Fig. 10(a), if $\\{\Lambda(\beta,\bar{f}),0\\}$, which is an
analytically obtained initial guess, resides in the region
$\Xi_{v}(v,v_{a})<0$, it is the solution of the problem (103). Otherwise, the
solution can be searched for iteratively, as illustrated in Fig 10(b), by
alternately searching in $v_{a}$ and $v$-directions from the initial guess
$\\{\Lambda(\beta,\bar{f}),0\\}$. The following is an algorithm to obtain the
solution:

$\displaystyle\mbox{Function }\Lambda_{\mathrm{reg}}(\beta,\bar{f})$ (105a)
$\displaystyle\quad v:=\Lambda(\beta,\bar{f})$ (105b)
$\displaystyle\quad\mbox{If
}u_{a}>0\,\wedge\,\hat{{\boldsymbol{u}}}\in{\boldsymbol{\mathcal{U}}}_{+}\,\wedge\,v>0\,\wedge\,\Xi_{v}(v,0)<0\
\mbox{Then }$ (105c) $\displaystyle\quad\quad v_{a}:=0$ (105d) Loop (105e)
$\displaystyle\quad\quad\quad
v_{a}:=\mbox{FindRoot}\left(\Xi_{v}(v,\bullet),[v_{a},v]\right)$ (105f)
$\displaystyle\quad\quad\quad\mbox{If }|\Xi_{f}(v,v_{a})|<\varepsilon_{f}\
\mbox{Exit Loop}$ (105g) $\displaystyle\quad\quad\quad
v:=\mbox{FindRoot}\left(\Xi_{f}(\bullet,v_{a}),[v,(F_{hM}-\bar{f})/\beta]\right)$
(105h) $\displaystyle\quad\quad\quad\mbox{If
}|\Xi_{v}(v,v_{a})|<\varepsilon_{v}\ \mbox{Exit Loop}$ (105i) While (105j) End
If (105k) $\displaystyle\mbox{Return }v.$ (105l)

Here, $\varepsilon_{f}$ and $\varepsilon_{v}$ are very small positive numbers.
This algorithm can be seen as an extension of $\Lambda(\beta,\bar{f})$ given
in Section 4.2. As is the case with $\Gamma_{\mathrm{reg}}$, the computation
of “FindRoot” can be performed with common iterative solvers. Some analytical
methods may be found because $\Xi_{f}$ is also a piece-wise parabolic
function.

Figure 9: Solution of the problem (103) in (a) the case of
$\Xi_{v}(\Lambda(\beta,\bar{f}),0)\geq 0$, where
$\\{v,v_{a}\\}=\\{\Lambda(\beta,\bar{f}),0\\}$ is the solution, and (b) the
case of $\Xi_{v}(\Lambda(\beta,\bar{f}),0)<0$, where iterative computation is
needed.

Figure 10: Simulation results of the system of Fig. 6 with the regeneration
valve being open ($u_{a}=0.5$) and closed ($u_{a}=0$).

Some simulations were performed with the same arm system as in Section 4.3 and
Fig. 6. The arm was moved down from $\theta=30$ deg for $1.5$ s with the
command $u_{c}=0.5$ with the regeneration valve being open ($u_{a}=0.5$) and
closed ($u_{a}=0$). The results are shown in Fig. 10. It shows that $u_{a}>0$
realizes a faster extending motion, which is consistent with the effects of
actual regeneration circuits. In addition, it can be seen that the impulsive
actuator force at the beginning of the motion (around $t=0.5$ s) is smaller
with $u_{a}>0$. It can be explained by the fact that the open regeneration
valve allows more flow of the oil, resulting in a softer response. The
impulses at the end of the movement (around $t=2$ s) should not be directly
compared because the angle and the velocity at this time are quite different
between the two conditions.

##  6 Extension 2: Multiple Actuators Driven by One Pump

###  6.1 Quasistatic Model

In some excavators, more than one actuators are actuated by a single pump as
illustrated in Fig. 11. In such a circuit, the behaviors of the actuators
influence each other, e.g., the movement of one actuator may decrease the
supplied flowrate to other actuators. This section presents an extension of
the quasistatic model to deal with such systems.

In circuits like Fig. 11, the actuators are connected to a junction at which
the pressure is $P$, and the total supplied flowrate $Q$ from the pump is
equal to the sum of the supplied flowrates to all actuators plus that
discharged through the bleed valve and the pump relief valve. On the other
hand, the actuator model developed in Section 3.3 assumes that the supplied
flowrate $Q$ from the pump is the input to be given. Therefore, one can see
that it is convenient to have a modified version of the quasistatic map
$\Gamma$ of which the input is the pressure $P$ instead of the flowrate $Q$.

Figure 11: Multiple actuators driven by one pump.

With a close look at the definition (64) of the function $\Gamma$ in Section
3.3, one can see that a $P$-input version of the quasistatic map $\Gamma$ can
be obtained by replacing $P_{m}$ by $P$ and removing the segments depending on
$Q$. Specifically, the modified version $\widehat{\Gamma}(P,v)$, which maps
the velocity $v$ to the force $f$ of an actuator depending on $P$, can be
given as follows:

$\displaystyle f$ $\displaystyle\in$
$\displaystyle\widehat{\Gamma}(P,v)\stackrel{{\scriptstyle\Delta}}{{=}}\mathrm{gsgn}(\widehat{\Gamma}_{-}(P,v),v,\widehat{\Gamma}_{+}(P,v))$
(106a) where $\displaystyle\widehat{\Gamma}_{+}(P,v)$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle\max(\min(\max(\Gamma_{+0a}(v),\Gamma_{+0b}(v)),$ (106b)
$\displaystyle\quad\max(\widehat{\Gamma}_{+2a}(P,v),\widehat{\Gamma}_{+2b}(P,v))),\Gamma_{+3}(v),-F_{rM})\quad$
$\displaystyle\widehat{\Gamma}_{-}(P,v)$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle\min(\max(\min(\Gamma_{-0a}(v),\Gamma_{-0b}(v)),$ (106c)
$\displaystyle\quad\min(\widehat{\Gamma}_{-2a}(P,v),\widehat{\Gamma}_{-2b}(P,v))),\Gamma_{-3}(v),F_{hM})\quad$
$\displaystyle\widehat{\Gamma}_{+2a}(P,v)$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$ $\displaystyle
A_{h}P-\dfrac{{\mathcal{S}}(v)}{\hat{u}_{ph}^{2}}-\dfrac{{\mathcal{S}}(v)}{\hat{u}_{tr}^{2}}$
(106d) $\displaystyle\widehat{\Gamma}_{+2b}(P,v)$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$ $\displaystyle
A_{h}P-\dfrac{{\mathcal{S}}(v)}{\hat{u}_{ph}^{2}}-F_{rM}$ (106e)
$\displaystyle\widehat{\Gamma}_{-2a}(P,v)$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$ $\displaystyle-
A_{r}P-\dfrac{{\mathcal{S}}(v)}{\hat{u}_{pr}^{2}}-\dfrac{{\mathcal{S}}(v)}{\hat{u}_{th}^{2}}$
(106f) $\displaystyle\widehat{\Gamma}_{-2b}(P,v)$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$ $\displaystyle-
A_{r}P-\dfrac{{\mathcal{S}}(v)}{\hat{u}_{pr}^{2}}+F_{hM}.$ (106g)

Here, the functions without hats are those defined in (64). In the same manner
as (65), we also have the following:

$\displaystyle\widehat{\Gamma}_{+}(P,0)=\widehat{\Gamma}_{h+}(P,0)-\widehat{\Gamma}_{r+}(P,0),\quad\widehat{\Gamma}_{-}(P,0)=\widehat{\Gamma}_{h-}(P,0)-\widehat{\Gamma}_{r-}(P,0)$
(107a) where $\displaystyle\widehat{\Gamma}_{h+}(P,0)$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}\min(F_{hM},A_{h}P)&\mbox{if
}\hat{u}_{ph}>0\\\ 0&\mbox{if }\hat{u}_{ph}=0\\\ \end{array}\right.$ (107d)
$\displaystyle\widehat{\Gamma}_{r+}(P,0)$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}0&\mbox{if }\hat{u}_{tr}>0\\\
F_{rM}&\mbox{if }\hat{u}_{tr}=0\\\ \end{array}\right.$ (107g)
$\displaystyle\widehat{\Gamma}_{h-}(P,0)$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}0&\mbox{if }\hat{u}_{th}>0\\\
F_{hM}&\mbox{if }\hat{u}_{th}=0\end{array}\right.$ (107j)
$\displaystyle\widehat{\Gamma}_{r-}(P,0)$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}\min(F_{rM},A_{r}P_{M})&\mbox{if
}\hat{u}_{pr}>0\\\ 0&\mbox{if }\hat{u}_{pr}=0.\\\ \end{array}\right.$ (107m)

The flowrate $Q_{p}$ into an actuator is determined by the pressure $P$ at the
junction and the pressure $P_{h}$ or $P_{r}$ of the chamber connected to the
pump, specifically, as follows:

$\displaystyle
Q_{p}=A_{h}\hat{u}_{ph}\max({\mathcal{R}}(A_{h}P-F_{h}),0)+A_{r}\hat{u}_{pr}\max({\mathcal{R}}(A_{r}P-F_{r}),0).$
(108)

If $\hat{{\boldsymbol{u}}}\in{\boldsymbol{\mathcal{U}}}_{+}$, we have
$F_{h}=f+F_{r}$ and
$F_{r}=\mathrm{proj}_{[0,F_{rM}]}({\mathcal{S}}(v)/\hat{u}_{tr}^{2})$. If
$\hat{{\boldsymbol{u}}}\in{\boldsymbol{\mathcal{U}}}_{-}$, we have
$F_{h}=\mathrm{proj}_{[0,F_{hM}]}({\mathcal{S}}(v)/\hat{u}_{th}^{2})$ and
$F_{r}=F_{h}-f$. Therefore, $Q_{p}$ is obtained as follows:

$\displaystyle Q_{p}=\widehat{Q}_{p}(P,v,\widehat{\Gamma}(P,v))$ (109)

where

$\displaystyle\widehat{Q}_{p}(P,v,\widehat{\Gamma}(P,v))$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle\left\\{\begin{array}[]{ll}A_{h}\hat{u}_{ph}\max({\mathcal{R}}(A_{h}P-\mathrm{proj}_{[0,F_{rM}]}({\mathcal{S}}(v)/\hat{u}_{tr}^{2})-f),0)&\mbox{if
}\hat{{\boldsymbol{u}}}\in{\boldsymbol{\mathcal{U}}}_{+}\\\ 0&\mbox{if
}\hat{{\boldsymbol{u}}}\in{\boldsymbol{\mathcal{U}}}_{0}\\\
A_{r}\hat{u}_{pr}\max({\mathcal{R}}(A_{r}P-\mathrm{proj}_{[0,F_{hM}]}(-{\mathcal{S}}(v)/\hat{u}_{th}^{2})+f),0)&\mbox{if
}\hat{{\boldsymbol{u}}}\in{\boldsymbol{\mathcal{U}}}_{-}.\\\
\end{array}\right.$ (113)

The function $\widehat{\Gamma}$ appearing in (109) may be set-valued at $v=0$,
but a careful observation of its limits of both sides of zero shows that
$\widehat{Q}_{p}(P,0,\widehat{\Gamma}(P,0))=0$, single-valued, under all three
conditions, ${\boldsymbol{\mathcal{U}}}_{+}$, ${\boldsymbol{\mathcal{U}}}_{0}$
and ${\boldsymbol{\mathcal{U}}}_{-}$.

By using the function $\widehat{Q}_{p}$, the algebraic constraint between the
total flowrate $Q$ from the pump and the pressure $P$ can be described as
follows:

$\displaystyle\Xi_{P}(P)\in{\mathcal{N}}_{(-\infty,P_{M}]}(P)$ (114a) where
$\displaystyle\Xi_{P}(P)$ $\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle
Q-U_{b}{\mathcal{R}}(P)-\sum_{j=1}^{N}\widehat{Q}_{p,j}(P,v_{j},\widehat{\Gamma}_{j}(P,v_{j})).$
(114b)

Here, the symbols with the subscript $j$ stand for those associated with the
$j$th actuator, and $N$ denotes the number of actuators. The value of
$\Xi_{P}(P)$ can be interpreted as the flowrate from the pump relief valve
(see Fig. 11), which is the difference between the oil supply from the pump
and the sum of the oil supplies to all actuators plus that to the bleed valve.
As long as $P<P_{M}$, (114a) reduces to $\Xi_{P}(P)=0$, which means that the
pump relief valve is closed. When $P=P_{M}$, (114a) reduces to $\Xi_{P}(P)\geq
0$, which means that the oil is discharged from the pump relief valve, of
which the pressure limit is $P_{M}$. The pressure $P$ can be found by solving
the algebraic problem (114) as follows:

$\displaystyle P$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}P_{M}&\mbox{if }\Xi_{P}(P_{M})\geq
0\\\
\mbox{FindRoot}(\Xi_{P}(\bullet),[0,P_{M}])&\mbox{otherwise.}\end{array}\right.$
(117)

This “FindRoot” is also easy because of the monotonicity of the function
$\Xi_{P}$.

By using the value of $P$ obtained by (117), the quasistatic relation between
the forces ${\boldsymbol{f}}=[f_{1},\cdots,f_{N}]^{T}$ and the velocities
${\boldsymbol{v}}=[v_{1},\cdots,v_{N}]^{T}$ of $N$ actuators are written in
the following form:

$\displaystyle{\boldsymbol{f}}\in{\boldsymbol{\Gamma}}_{\mathrm{mul}}({\boldsymbol{v}})\stackrel{{\scriptstyle\Delta}}{{=}}\left[\widehat{\Gamma}_{1}(P,v_{1}),\cdots,\widehat{\Gamma}_{N}(P,v_{N})\right]^{T}$
(118)

where $P$ is the one obtained by (117).

Some numerical examples are shown in Fig. 12. In these examples, two identical
actuators 1 and 2 share a single pump. The parameters of the actuators are the
same as those in Section 3.4. The force $f_{1}$ obtained by the map
$[f_{1},f_{2}]^{T}\in{\boldsymbol{\Gamma}}_{\mathrm{mul}}([v_{1},v_{2}]^{T})$
according to the variable $v_{1}$ and some fixed values of $v_{2}$ are shown
in Figs. 12(a) and (d). Intermediate values $P$ and $\Xi_{P}(P)$, which are
the immediate output of the root finding in (117), are presented in Figs.
12(b)(c) and (e)(f). It can be seen that an increased speed $v_{2}$ of the
actuator 2 results in a decreased speed $v_{1}$ of the actuator 1, a decreased
pressure $P$ at the junction, and a decreased flowrate $\Xi_{P}(P)$ from the
pump relief valve, which are consistent with what can happen in real hydraulic
circuits.

Figure 12: Numerical examples regarding the quasistatic map
$[f_{1},f_{2}]^{T}\in{\boldsymbol{\Gamma}}_{\mathrm{mul}}([v_{1},v_{2}]^{T})$
of a circuit including two identical actuators driven by a single pump.
Positive commands $u_{c,1}=u_{c,2}=0.5$ are given in (a)-(c) and negative
commands $u_{c,1}=u_{c,2}=-0.5$ are given in (d)-(f) (see (67)).

###  6.2 Incorporation into Simulators

In the same way as in Section 4 and in Section 5.2, the application of the
quasistatic map ${\boldsymbol{\Gamma}}_{\mathrm{mul}}$ to multibody
simulations requires a map ${{\boldsymbol{\Lambda}}}_{\mathrm{mul}}$ that
satisfies the following relation:

$\displaystyle\mathrm{diag}{[}{\boldsymbol{\beta}}]{\boldsymbol{v}}+\bar{{\boldsymbol{f}}}\in{\boldsymbol{\Gamma}}_{\mathrm{mul}}({\boldsymbol{v}})$
$\displaystyle\iff$
$\displaystyle{\boldsymbol{v}}\in{{\boldsymbol{\Lambda}}}_{\mathrm{mul}}({\boldsymbol{\beta}},\bar{{\boldsymbol{f}}})$
(119)

where ${\boldsymbol{\beta}}$, ${\boldsymbol{v}}$ and $\bar{{\boldsymbol{f}}}$
are $N$-dimenional vectors and the elements of ${\boldsymbol{\beta}}$ are all
positive. The structure (118) of ${\boldsymbol{\Gamma}}_{\mathrm{mul}}$
suggests that ${{\boldsymbol{\Lambda}}}_{\mathrm{mul}}$ has the following
structure:

$\displaystyle{{\boldsymbol{\Lambda}}}_{\mathrm{mul}}({\boldsymbol{\beta}},\bar{{\boldsymbol{f}}})=[\widehat{\Lambda}_{1}(P,\beta_{1},\bar{f}_{1}),\cdots,\widehat{\Lambda}_{N}(P,\beta_{N},\bar{f}_{N})]^{T}$
(120)

where $\widehat{\Lambda}_{j}(P,\beta_{j},\bar{f}_{j})$ are those that satisfy
the following:

$\displaystyle\beta_{j}v_{j}+\bar{f}_{j}\in\widehat{\Gamma}_{j}(P,v_{j})\iff
v_{j}\in\widehat{\Lambda}(P,\beta_{j},\bar{f}).$ (121)

Recalling that $\Lambda$ is obtained as (73) and (74), one can obtain
$\widehat{\Lambda}_{j}$ as follows:

$\displaystyle\widehat{\Lambda}(P,\beta,\bar{f})$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}v_{hM}&\mbox{if }\bar{f}\geq
F_{hM}\wedge{\boldsymbol{u}}\in{\boldsymbol{\mathcal{U}}}_{+}\\\
\max(\min(\max(v_{+0a},v_{+0b}),\max(\widehat{v}_{+2a},\widehat{v}_{+2b})),v_{+3},v_{rM})&\\\
&\mbox{if
}\bar{f}\leq\widehat{\Gamma}_{+}(P,0)\wedge{\boldsymbol{u}}\in{\boldsymbol{\mathcal{U}}}_{+}\\\
v_{rM}&\mbox{if }\bar{f}\leq-F_{rM}\ \wedge\
{\boldsymbol{u}}\in{\boldsymbol{\mathcal{U}}}_{-}\\\
\min(\max(\min(v_{-0a},v_{-0b}),\min(\widehat{v}_{-2a},\widehat{v}_{-2b})),v_{-3},v_{hM})&\\\
&\mbox{if }\bar{f}\geq\widehat{\Gamma}_{-}(P,0)\ \wedge\
{\boldsymbol{u}}\in{\boldsymbol{\mathcal{U}}}_{-}\\\ 0&\mbox{if
}\widehat{\Gamma}_{+}(P,0)<\bar{f}<\widehat{\Gamma}_{-}(P,0)\end{array}\right.$
(122h) where $\displaystyle\widehat{v}_{+2a}$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle\Phi_{A}\left(\beta,\bar{f}-A_{h}P,\psi(\hat{u}_{ph},\hat{u}_{tr})\right)$
(122i) $\displaystyle\widehat{v}_{+2b}$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle\Phi_{A}\left(\beta,\bar{f}-A_{h}P+F_{rM},\hat{u}_{ph}^{2}\right)$
(122j) $\displaystyle\widehat{v}_{-2a}$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle\Phi_{A}\left(\beta,\bar{f}+A_{r}P,\psi(\hat{u}_{pr},\hat{u}_{th})\right)$
(122k) $\displaystyle\widehat{v}_{-2b}$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$
$\displaystyle\Phi_{A}\left(\beta,\bar{f}+A_{r}P-F_{hM},\hat{u}_{pr}^{2}\right).$
(122l)

Here, the subscript $j$ denoting the actuator indices are omitted and the
symbols without hats are those defined in (73).

Figure 13: Simulation results of two arms identical to the one in Fig. 6
driven by a single pump. No external forces were applied to either arm. The
command $u_{c,1}$ sent to the arm 1 was the same between the two scenarios,
while the arm 2 was actuated only in Scenario II.

In the expression (120), the pressure $P$ can be obtained in the same line of
thought as in Section 6.1, but the problem (114) needs to be modified because
$v_{j}$ are unknowns. Considering that $\bar{{\boldsymbol{f}}}$ is given, the
problem is redefined as

$\displaystyle\Xi_{P\Lambda}(P)\in{\mathcal{N}}_{(-\infty,P_{M}]}(P)$ (123a)
where $\displaystyle\Xi_{P\Lambda}(P)$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}$ $\displaystyle
Q-U_{b}{\mathcal{R}}(P)-\sum_{j=1}^{N}\widehat{Q}_{p,j}(P,\widehat{\Lambda}_{j}(\beta_{j},\bar{f}_{j}),\bar{f}_{j}+\beta_{j}\widehat{\Lambda}_{j}(\beta_{j},\bar{f}_{j})).$
(123b)

The solution $P$ of the problem (123) can be written as follows:

$\displaystyle P$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}P_{M}&\mbox{if
}\Xi_{P\Lambda}(P_{M})\geq 0\\\
\mbox{FindRoot}(\Xi_{P\Lambda}(\bullet),[0,P_{M}])&\mbox{otherwise.}\end{array}\right.$
(126)

In conclusion, the function ${{\boldsymbol{\Lambda}}}_{\mathrm{mul}}$
satisfying (119) is obtained as (120) in which $\widehat{\Lambda}_{j}$ are
those defined by (122) and $P$ is replaced by (126).

Some simulations were performed with a system composed of two arms identical
to the one in Fig. 6 connected with a single pump. It was implemented with the
function ${{\boldsymbol{\Lambda}}}_{\mathrm{mul}}$ with $N=2$. The parameters
of the arms were set the same as in Section 4.3, and the parameters of the
actuators and the circuit were set the same as in Section 6.1. No external
forces were applied to either arm. The results are shown in Fig. 13. The
initial postures of the arms 1 and 2 were set as $\theta=-30$ deg and
$\theta=20$ deg, respectively. In the first scenario, the arm 1 was driven
upward by the command $u_{c,1}=-0.6$ from $t=0.5$ s to $2.5$ s. while the arm
2 was not driven. In the second scenario, the command to the arm 1 was the
same as that in the first scenario, the arm 2 was driven downward from $t=1$ s
to $2$ s. As can be seen in Figs. 13(b) and (c), the actuation of the arm 2
affected the motion of the arm 1, i.e., the arm 1 was decelerated when the arm
2 was actuated. It can be explained by the effect of the arm 2, which
suctioned the oil when actuated, causing the shortage of the oil flow into the
arm 1.

##  7 Conclusions

This article has presented a quasistatic model of a hydraulic actuator driven
by a four-valve independent metering circuit. The presented model is described
as a nonsmooth map between the velocity and the force. The model is derived
from the algebraic constraint between the flowrate and the pressure at every
valve in the circuit in the steady state. This article also presents an
approach to incorporate the quasistatic model into the multibody dynamics
simulators, in which the hydraulic model is connected to rigid bodies through
virtual viscoelastic elements. In addition, the proposed model is extended to
include a regeneration pipeline and to deal with a collection of actuators
driven by a single pump.

In multibody simulations employing the presented quasistatic model, the
transient responses are determined by the virtual viscoelastic elements, which
are parameterized by the stiffness $K$ and the viscosity $B$. In reality,
transient responses are governed by many factors such as the compressibility
of the oil, the inertia of the actuator and the oil, and the compliance of the
pipes. Some tuning guidelines for the parameters $K$ and $B$ should be sought
to reproduce the behaviors of actual systems by the presented simulation
framework. In order to accelerate the computation, analytical methods for the
root-finding routines that have appeared in the proposed algorithm in Sections
5 and 6 should also be addressed. Integration of the schemes of Sections 5 and
6, i.e., multiple actuators with regeneration pipelines driven by a single
pump, is also an open problem.

## References

  * [1] V. Acary, O. Bonnefon, and B. Brogliato.  Time-stepping numerical simulation of switched circuits within the nonsmooth dynamical systems approach.  IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, 29(7):1042–1055, 2010.
  * [2] V. Acary and B. Brogliato.  Numerical Methods for Nonsmooth Dynamical Systems: Applications in Mechanics and Electronics, volume 35 of Lecture Notes in Applied and Computational Mechanics.  Springer, 2008.
  * [3] K. Addi, S. Adly, B. Brogliato, and D. Goeleven.  A method using the approach of Moreau and Panagiotopoulos for the mathematical formulation of non-regular circuits in electronics.  Nonlinear Analysis: Hybrid Systems, 1(1):30–43, 2007.
  * [4] W. Borutzky, B. Barnard, and J. Thoma.  An orifice flow model for laminar and turbulent conditions.  Simulation Modelling Practice and Theory, 10:141–152, 2002.
  * [5] B. Brogliato, A. Daniilidis, C. Lemaréchal, and V. Acary.  On the equivalence between complementarity systems, projected systems and differential inclusions.  Systems & Control Letters, 55(1):45–51, 2006.
  * [6] K. Choi, J. Seo, Y. Nam, and K. U. Kim.  Energy-saving in excavators with application of independent metering valve.  Journal of Mechanical Science and Technology, 29(1):387–395, 2015\.
  * [7] D. Christofori and A. Vacca.  Modeling hydraulic actuator mechanical dynamics from pressure measured at control valve ports.  Proceedings of the Institution of Mechanical Engineers, Part I: Journal of Systems and Control Engineering, 229(6):541–558, 2015.
  * [8] M. C. Destro and V. J. De Negri.  Method for combining valves with symmetric and asymmetric cylinders for hydraulic systems.  International Journal of Fluid Power, 19(3):126–139, 2018.
  * [9] B. Eriksson and J.-O. Palmberg.  Individual metering fluid power systems: challenges and opportunities.  Proceedings of IMechE: Journal of Systems and Control Engineering, 225(2):196–211, 2011.
  * [10] A. Esposito.  A simplified method for analyzing hydraulic circuits by analogy.  Machine Design, 41(24):173–177, 1969.
  * [11] D. Goeleven.  Existence and uniqueness for a linear mixed variational inequality arising in electrical circuits with transistors.  Journal of Optimization Theory and Applications, 138:397–406, 2008.
  * [12] J. Greenslade, Thomas B.  The hydraulic analogy for electric current.  Physics Teacher, 41:464–466, 2003.
  * [13] M. Kiani Oshtorjani, A. Mikkola, and P. Jalali.  Numerical treatment of singularity in hydraulic circuits using singular perturbation theory.  IEEE/ASME Transactions on Mechatronics, 24(1):144–153, 2018.
  * [14] R. Kikuuwe.  A sliding-mode-like position controller for admittance control with bounded actuator force.  IEEE/ASME Transactions on Mechatronics, 19(5):1489–1500, 2014\.
  * [15] R. Kikuuwe.  Anti-noise and anti-disturbance properties of differential-algebraic relaxation applied to a set-valued controller.  In Proceedings of the 15th International Workshop on Variable Structure Systems (VSS18), pages 480–485, 2018.
  * [16] R. Kikuuwe.  A brush-type tire model with nonsmooth representation.  Mathematical Problems in Engineering, 2019:9747605, 2019.
  * [17] R. Kikuuwe, N. Takesue, A. Sano, H. Mochiyama, and H. Fujimoto.  Admittance and impedance representations of friction based on implicit Euler integration.  IEEE Transactions on Robotics, 22(6):1176–1188, 2006.
  * [18] R. Kikuuwe, S. Yasukouchi, H. Fujimoto, and M. Yamamoto.  Proxy-based sliding mode control: A safer extension of PID position control.  IEEE Transactions on Robotics, 26(4):670–683, 2010.
  * [19] A. Lichtarowicz, R. K. Duggins, and E. Markland.  Discharge coefficients for incompressible non-cavitating flow through long orifices.  Journal of Mechanical Engineering Science, 7(2):210–219, 1965\.
  * [20] J. Rahikainen, F. González, M. Á. Naya, J. Sopanen, and A. Mikkola.  On the cosimulation of multibody systems and hydraulic dynamics.  Multibody System Dynamics, 50(2):143–167, 2020.
  * [21] J. Rahikainen, M. Kiani, J. Sopanen, P. Jalali, and A. Mikkola.  Computationally efficient approach for simulation of multibody and hydraulic dynamics.  Mechanism and Machine Theory, 130:435–446, 2018.
  * [22] J. Rahikainen, A. Mikkola, J. Sopanen, and J. Gerstmayr.  Combined semi-recursive formulation and lumped fluid method for monolithic simulation of multibody and hydraulic dynamics.  Multibody System Dynamics, 44:293–311, 2018.
  * [23] M. Ruderman.  Full- and reduced-order model of hydraulic cylinder for motion control.  In Proceedings of 43rd Annual Conference of IEEE Industrial Electronics Society (IECON 2017), pages 7275–7280, 2017.
  * [24] S. Sakai and S. Stramigioli.  Visualization of hydraulic cylinder dynamics by a structure preserving nondimensionalization.  IEEE/ASME Transactions on Mechatronics, 23(5):2196–2206, 2018\.
  * [25] A. Shenouda.  Quasi-static hydraulic control systems and energy savings potential using independent metering four-valve assembly configuration.  PhD thesis, Georgia Institute of Technology, 2006.
  * [26] A. Shenouda and W. Book.  Optimal mode switching for a hydraulic actuator controlled with four-valve independent metering configuration.  International Journal of Fluid Power, 9(1):35–43, 2008.
  * [27] K. A. Tabor.  A novel method of controlling a hydraulic actuator with four valve independent metering using load feedback.  SAE Technical Paper Series, 2005-01-3639, 2005.
  * [28] L. Wang, W. J. Book, and J. D. Huggins.  Application of singular perturbation theory to hydraulic pump controlled systems.  IEEE/ASME Transactions on Mechatronics, 17(2):251–259, 2012.
  * [29] D. Wu, R. Burton, and G. Schoenau.  An empirical discharge coefficient model for orifice flow.  International Journal of Fluid Power, 3(3):13–10, 2002.
  * [30] X. Xiong, R. Kikuuwe, and M. Yamamoto.  A differential-algebraic method to approximate nonsmooth mechanical systems by ordinary differential equations.  Journal of Applied Mathematics, 2013:320276, 2013.
  * [31] X. Xiong, R. Kikuuwe, and M. Yamamoto.  A multistate friction model described by continuous differential equations.  Tribology Letters, 51(3):513–523, 2013.
  * [32] J. Yao, Z. Jiao, D. Ma, and L. Yan.  High-accuracy tracking control of hydraulic rotary actuators with modeling uncertainties.  IEEE/ASME Transactions on Mechatronics, 19(2):633–641, 2014.
  * [33] Y. Ye, C.-B. Yin, X.-D. Li, W.-J. Zhou, and F.-F. Yuan.  Effects of groove shape of notch on the flow characteristics of spool valve.  Energy Conversion and Management, 86:1091–1101, 2014.
  * [34] A. Ylinen, H. Marjamäki, and J. Mäkinen.  A hydraulic cylinder model for multibody simulations.  Computers and Structures, 138(1):67–72, 2014.