WGMwgm QEqe EPep PMSpms BECbec DEde

[1] [1] [2]

[cor2]These Authors equally contributed to the Manuscript [cor1]Corresponding
author

# A compartmental model for cyber-epidemics

D. Aleja  G. Contreras-Aso  K. Alfaro-Bittner  E. Primo<EMAIL_ADDRESS>R.
Criado  M. Romance  S. Boccaletti  Universidad Rey Juan Carlos, Calle Tulipán
s/n, Móstoles, 28933, Madrid,Spain  Laboratory of Mathematical Computation on
Complex Networks and their Applications (LaCoMaRCA), Universidad Rey Juan
Carlos, Calle Tulipán s/n, Móstoles, 28933 , Madrid, Spain  Data, Complex
networks and Cybersecurity Research Institute, Universidad Rey Juan Carlos,
Plaza Manuel Becerra 14, 28028, Madrid, Spain  Moscow Institute of Physics and
Technology, Dolgoprudny, 141701, Moscow, Russian Federation  CNR - Institute
of Complex Systems, Via Madonna del Piano 10, I-50019, Sesto Fiorentino, Italy

###### Abstract

In our more and more interconnected world, a specific risk is that of a cyber-
epidemic (or cyber-pandemic), produced either accidentally or intentionally,
where a cyber virus propagates from device to device up to undermining the
global Internet system with devastating consequences in terms of economic
costs and societal harms related to the shutdown of essential services. We
introduce a compartmental model for studying the spreading of a malware and of
the awareness of its incidence through different waves which are evolving on
top of the same graph structure (the global network of connected devices).
This is realized by considering vectorial compartments made of two components,
the first being descriptive of the state of the device with respect to the new
malware’s propagation, and the second accounting for the awareness of the
device’s user about the presence of the cyber threat. By introducing suitable
transition rates between such compartments, one can then follow the evolution
of a cyber-epidemic from the moment at which a new virus is seeded in the
network, up to when a given user realizes that his/her device has suffered a
damage and consequently starts a wave of awareness which eventually ends up
with the development of a proper antivirus software. We then compare the
overall damage that a malware is able to produce in Erdős-Rényi and scale-free
network architectures for both the case in which the virus is causing a fixed
damage on each device and the case where, instead, the virus is engineered to
mutate while replicating from device to device. Our result constitute actually
the attempt to build a specific compartmental model whose variables and
parameters are entirely customized for describing cyber-epidemics.

###### keywords:

Cyber-epidemics Compartmental model Networks Epidemics spreading

The study of compartmental models started already at the beginning of the 20th
century [1, 2, 3, 4], and soon became a subject of great, recent, interest
that has attracted the attention of many epidemiologists. In analogy with
cellular automata [5], these models consider a networked population of
individuals, each one of them described by a state whose discrete values are
labeled by compartments. Individuals may then progress between compartments
through given transition probabilities, which allow the time-discrete
evolution of the population during, for instance, the spreading of infectious
diseases and/or rumors and social contagion. Physicists became interested in
these models when it was pointed out that epidemiological processes can be
regarded as percolation like processes [6].

Starting from the seminal work by Pastor–Satorras and Vespignani [7], the last
twenty years have seen a burst of activity on understanding the effects of a
network topology on the rate and patterns of the disease spread. A lot of
studies have tried indeed to predict things such as the total number of
infected individuals, or the duration of an epidemic, and to estimate various
relevant parameters such as the reproductive number. Moreover, and especially
in relation to the recent COVID19 world pandemic crisis, these models have
been used to assess the effects of public health interventions and/or to
quantify the efficiency of issuing a limited number of vaccines in a given
population [8, 9].

In particular, the SIR model [10, 11] is one of the simplest compartmental
models, and describes the evolution of diseases resulting in the immunization
or death of the infected individuals. The model assumes that, at each time,
each individual can be in one of three possible compartments: susceptible
(denoted by $S$), infected ($I$), or removed ($R$). The susceptible units of
the network are those healthy persons that can develop the disease if they get
in contact with infected individuals. Once an individual contracts the
infection, it moves into the infected (and infective) compartment, and then,
after some time, into the removed compartment, which indicates that the
individual cannot catch the disease anymore (or passes it on), due to a
lasting resistance conferred by the recovery (or because it dies).

Figure 1: The compartmental model for cyber-epidemics. Units of the network
can be in one of four possible states before recovery: susceptible-unaware
($S_{u}$), susceptible-aware ($S_{a}$), infected-unaware ($I_{u}$), and
infected aware ($I_{a}$). The recovery or healed status ($H_{a}$) is by nodes
that, after getting aware of the existence of the malware, applies a suitable
operation on the system. Orange (blue) arrows stand for contact-based
(localized spontaneous) transitions.

Just as a viral pandemic proliferates inside a population of individuals, a
ransomware (or other malware) software can spread within the global Internet.
As technology is today globally interconnected, a cyber virus can indeed
propagate from device to device, with devastating consequences. Cyber attacks
on infrastructure services are currently on the rise, and hackers are
exploiting the use of the Internet of Things which creates millions of new
vulnerability points in all critical infrastructures. A specific risk is
therefore the occurrence of a cyber-epidemics (or cyber-pandemic), produced
either accidentally or intentionally, and undermining the global Internet
system up to the need of its lockdown. The World Economic Forum [12] predicted
that a single day without Internet would cost around 50 billion USD globally,
without even considering the societal harm related to the shutdown of
essential services.

Since the first attempt to describe the spreading of computer viruses [13],
several other studies have tried to adapt classical, global or networked,
compartmental models with the aim not only of investigating the propagation
cycle of cyber viruses but also to evaluate the effectiveness of possible
security countermeasures [14, 15, 16, 17, 18, 19, 20]. There is, however, a
fundamental difference between the spreading of a biological virus within a
population and a cyber-epidemic. In the former case, indeed, each individual
is passive actor of the game, contracting the disease and recovering from it
due to the action of its immune system. The latter case can instead be seen as
the struggle between two kind of actors: the ones who intentionally program
the malicious code and try to seed it within the global Internet (that we will
call from here after as the “bad team”) and those (that we will call from here
after as the “good team”) who are instead engaged in programming the
corresponding antivirus code after becoming aware of the presence of the new
malware in the network, and in spreading it to all network’s users.

In this paper, we introduce a novel compartmental model, whose variables and
parameters are entirely customized for the case of a cyber-epidemic. We will
then compare the spreading of a malware on top of Erdős-Rényi (connected) and
scale-free network architectures for both the case in which the virus is
producing a fixed damage on each device and the case where, instead, the virus
is engineered to mutate while replicating from device to device.

In our model, both the spreading of malware and awareness occurs through waves
developing and evolving within the same graph structure (the global network of
connected devices). Notice that this differentiates our approach from that of
Ref. [21] which describes the dynamical interplay of a virus and of an
awareness level on top of multiplex networks. In our case, instead, we
consider two dimensional vectorial compartments, the first component of which
being the state of the device with respect to the new malware’s propagation
($S$ for susceptible, $I$ for infected, and $H$ for healed) while the second
component accounts for the awareness of the device’s user about the presence
of the new virus (“a” for aware and “u” for unaware). Transition rates between
such vectorial compartments are then defined to properly model the behavior of
a cyber-epidemic from the moment at which a new (i.e., yet unknown to the good
team) virus is seeded in the network by the bad team, up to the moment at
which a given user realizes that his/her device has suffered a damage and
consequently starts a wave of awareness which eventually ends up with reaching
one of the units of the good team. Then, the good team develops a proper
antivirus software that the aware units download from the net for healing
their devices.

In other words, the framework we adopt is somewhat like a vectorial version of
the SIR model, where each susceptible or infected node is also either aware or
unaware of the existence of the virus. This means that one has four possible
states before recovery: susceptible-unaware ($S_{u}$), susceptible-aware
($S_{a}$), infected-unaware ($I_{u}$), and infected aware ($I_{a}$), see Fig.
1. Finally, the recovery or healed status ($H_{a}$) can be reached only by a
node (a device) whose user, after getting aware of the existence of the virus,
applies a suitable operation on the system, for instance installing an
antivirus.

The transition probabilities between such states are mediated by four
fundamental parameters, whose meaning is directly linked to specific processes
occurring during a cyber epidemic:

  * •

$\tau$ is the equivalent of the standard infection parameter in the SIR model.
Here, it accounts for the contact-based transition rate from the susceptible
to the infected status. In the case of $S_{a}$ individuals, a different
infection rate is used ($\tau^{\prime}=\tau/10$), due to the fact that
awareness makes individuals reluctant to be in contact with infected
sites/computers (for instance, email phishing campaigns normally raise alarms
in media, making people more concerned about clicking on dubious links);

  * •

$\nu$ is the rate of spreading of awareness due to network contacts, i.e., the
parameter who rules the contact-based transition from the unaware to the aware
status;

  * •

$\mu$ accounts for the individual awareness parameter, i.e., the rate at which
an infected unaware user actually notices that his/her device has suffered a
damage (larger than a given threshold) due to having been infected by a new
cyber virus, and consequently turns to the infected aware state and
simultaneously starts signaling the existence of a new virus circulating in
the net;

  * •

$\gamma$ is the recovery/healing parameter, i.e., the rate at which an aware
(either susceptible or infected) individual gets healed by purging of the
system via the installation of a proper antivirus.

Symbol | Description | Range | Chosen value
---|---|---|---
$\tau$ | Infection rate | [0,1] | 0.0055
$\nu$ | (Contact-based) awareness parameter | [0,1] | 0.011
$\mu_{0}$ | (Spontaneous or local) awareness parameter | [0,1] | 0.011
$\gamma$ | Recovery parameter | [0,1] | 0.03
$\rho_{0}$ | Fraction of population initially infected | [0,1] | 0.01
Table 1: Summary of the parameters involved in the model.

Two types of transitions may therefore occur in our model: contact-based
transitions and individual (or contact-independent) ones. In individual
transitions, the change of state is fully independent of the states of the
rest of the devices, but is only due to the perception of the damage caused to
the device (in the case of the passage from the state $I_{u}$ to the state
$I_{a}$ at rate $\mu$) and/or the user’s interest in using an antivirus (in
the case of the passage from any aware state to the state $H_{a}$ at rate
$\gamma$).

Figure 2: The typical cycle of a cyber-epidemic. The total number of $S_{u}$,
$S_{a}$, $I_{u}$, $I_{a}$ and $H_{a}$ units (see legend for color code) vs.
time, during a typical cyclic evolution of the virus. Simulations refer to an
ensemble average over 1000 realization of an Erdős-Rényi random network with
$N=1000$ nodes, average degree $\langle k\rangle=10$, threshold $\theta=0.2$
and constant damage $d=0.3$. See Table 1 for the values of the other
parameters used. It is possible to distinguish three different phases of the
cycle: an initial phase of no awareness where more and more nodes pass from
the $S_{u}$ to the $I_{u}$ state; a second stage, the awareness phase starting
at around $t\sim 30$, when a given $I_{u}$ node becomes $I_{a}$ and ignites an
awareness wave so that the total number of both $I_{a}$ and $S_{a}$ grows by
contact transitions; a final healing stage where all nodes turn to the $H_{a}$
state.

In our simulations, we initially prepare the network with all its units in the
$S_{u}$ state. Then, a new malware is seeded by the bad team, i.e., one node
(or a small group of nodes) is turned to the state $I_{u}$, and contact
propagations start, yielding an initial spreading of the virus. At a second
time, i.e., when one of the infected users gets aware of the damage produced
by the cyber virus in its device, a second wave (which spreads awareness in
the same network) starts due to an initial local transition from the state
$I_{u}$ to the state $I_{a}$. Eventually, the entire cycle of the cyber-
epidemic takes place, with an end state of the network where all its units are
in the $H_{a}$ state [22].

In order to properly monitor how damaging a cyber-epidemic can be in a network
of devices, we introduce as a further parameter the damage $d\in[0,1]$ caused
when a device is infected by the virus, and we quantify the total damage
caused to the system as the sum of the damages made to each single device.
Finally, as the $\mu$ parameter is totally related to the damage produced in a
device, it is appropriate to define it to be proportional to the $d$
parameter. In addition, in order to account for the fact that there may be
situations where the virus damage is not perceived at all by the user and the
device continues to function correctly even with high damage or, on the
opposite, cases (like, for instance, in high sensitivity infrastructures)
where just a little damage constitutes an alarm on the existence of the virus,
the parameter $\mu$ is activated depending on whether it exceeds a threshold
$\theta\in[0,1]$. The latter lead to the following expression with
$\mu_{0}\in[0,1]$:

$\mu=\begin{cases}\displaystyle\mu_{0}(d-\theta)\;\;\hbox{if}\;\;d\geq\theta,\\\
\displaystyle\mu=0\;\;\hbox{if}\;\;d<\theta.\end{cases}$

As for the threshold $\theta$, the higher it is the more damage individual
systems can withstand before realizing the presence of the virus. If the
threshold is at a low value (as, for instance, in high security
infrastructures or clusters) a damage as small as the deletion of a couple of
files in the system is already sufficient to trigger an alarm. If, instead,
the threshold is high, then a few files missing would be taken as an error,
but a large system failure would actually be understood as the consequence of
a virus. Lower thresholds are therefore always preferable for protecting an
overall system from damage, but in practice they are overcostly as they need
more investments in infrastructure protection and even more employees to
maintain it.

Fig. 2 reports the total number of $S_{u}$, $S_{a}$, $I_{u}$, $I_{a}$ and
$H_{a}$ units vs. time during a typical cyclic evolution of the virus, for an
Erdős-Rényi random network of $N=1000$ nodes, with average degree $\langle
k\rangle=10$, $\theta=0.2$, $\tau=0.0055$, $\nu=0.011$, $\mu_{0}=0.011$,
$\gamma=0.03$ and $\rho_{0}=0.01$ (being $\rho_{0}$ the fraction of population
initially infected, which therefore means that the virus is initially seeded
in 10 devices). From the figure, one can clearly distinguish the three
different stages of the cycle: an initial stage of no awareness where more and
more nodes are infected (passing from the $S_{u}$ to the $I_{u}$ state), the
beginning of the awareness phase (around $t\sim 30$) when a given $I_{u}$ node
becomes $I_{a}$ and starts an awareness wave so that the total number of both
$I_{a}$ and $S_{a}$ grows by contact transitions, and the final healing stage
where all nodes turn to the $H_{a}$ state.

In order to quantify the total injury produced in the system by the cyber
virus during its cycle of evolution, we introduce the quantity $D/N$
accounting for the normalized sum of the individual damages suffered by each
node. Its value is obtained with multiplying $d$ by the total number of
infected nodes during a cycle (regardless on whether they are in the state
$I_{u}$ or $I_{a}$), and dividing by $N$. When $d$ is below the threshold
$\theta$ the presence of the cyber threat is never detected and, as a
consequence, the virus will propagate to all nodes in the network yielding
$D/N=(N\cdot d)/N=d$. A non trivial behavior is instead observed for all
values $d>\theta$, where the awareness mechanism is activated at a given time
in the cycle.

As a first step, we face the problem of assessing how different topologies of
the network react to a virus (causing a constant damage $d$ when infecting a
device) at different values of the threshold $\theta$.

Figure 3: Effects of network topology. $D/N$ vs. $d$, for different choices of
networks (Barabási-Albert or Erdős-Rényi) of mean degree 6, 10 and 14 (see
legend at the bottom of the figure for color code) and different values of the
threshold $\theta$: a) $\theta=0.2$, b) $\theta=0.4$, and c) $\theta=0.6$.

With the aim of comparing homogeneous and heterogeneous topologies Fig. 3
illustrates the results of our simulations, and reports $D/N$ vs. $d$ for
Barabási-Albert [23, 24] scale free and Erdős-Rényi [25] random networks, at
different values of the mean degree $\langle k\rangle$ and different values of
$\theta$. A first, even though rather trivial, evidence is that the global
damage is higher for higher values of $\theta$, indicating that the harder it
is for nodes to become aware, the more free is the cyber threat to propagate
without countermeasures. Moreover, a noticeable difference is observed
regarding the position of the maximum damage that a virus can deal in a given
system, being it located at $d=1$ for low thresholds and at $d=\theta$ for
high thresholds. Finally, and more remarkably, the results shown in Fig. 3
allow to conclude that Barabási-Albert scale free (i.e., heterogenous)
networks are more fragile, and more sensitive to the spreading of viruses
causing a fixed damage than Erdős-Rényi (i.e., homogenous) ones, so that a
first conclusion is that engineering a network in a scale free topology
renders it more vulnerable to cyber-epidemics, in line with what was already
known about their structural fragility against intentional attacks [26]. On
the other hand, it is seen that, regardless on the specific network topology,
the higher is the average degree the higher is the graph’s vulnerability.
Already at this qualitative level, a first conclusion can be drawn: if the
threshold is high, inflicting maximal injury to a system requires some tuning
on the side of the bad team (the virus designers), with the ideal value of the
damage to be caused at each single device being that of the threshold. On the
opposite, if the threshold is low, then a high base damage virus is able to
inflict a huge damage on the network regardless on the specific threshold
value.

Figure 4: Effects of time varying damage. $D/N$ vs. $\varepsilon$ (see text
for definition), for different choices of networks (Erdős-Rényi in the left
column and Barabási-Albert in the right column) of mean degree 6, 10 and 14
(see legend at the bottom of the figure for color code) and different values
of the threshold $\theta$: a,b) $\theta=0.2$, c,d) $\theta=0.6$. In all panels
the maximal injuries achieved by a fixed strength virus (the global maximum of
the curves of Fig. 3) are reported with horizontal dashed lines.

We then move to consider the effects of a damage which is variable in time,
that is, $d=d(t)$, where $t$ is a discrete time measuring the number of
iterations in the model (a time unit being the lapse from one to another
iteration of the networks’ nodes). Namely, starting at $d_{0}=0.1$, the damage
is increased in times as

$d(t)=\frac{d_{0}e^{\varepsilon t}}{1+d_{0}(e^{\varepsilon t}-1)},$ (1)

for some $\varepsilon>0$. Note that the function $d(t)$ is the logistic
function with growth rate $\varepsilon$, initial population $d_{0}$ and
carrying capacity $1$, so that it displays an exponential behavior
($d(t)\approx d_{0}e^{\varepsilon t}$) when $t\approx 0$ and approaches $1$
when $t\to\infty$. Choosing the logistic function as the damage function is
ideal in the model, as the intentions of the bad team is precisely that of
spawning as much damage as possible before the virus is detected by users.

Fig. 4 reports $D/N$ vs. $\varepsilon$, for both Erdős-Rényi (left column) and
Barabási-Albert (right column) at different mean degrees and different values
of $\theta$. For comparison, in the same figure the maximum possible injury
caused to the system by a fixed strength virus is reported by horizontal
dashed lines. One can easily see that, as a function of $\varepsilon$, $D/N$
displays an initial monotonic growth and an asymptotic behavior for
$\varepsilon\to 1$, featuring a local maximum for some
$0<\bar{\varepsilon}<1$. Remarkably, one can notice that for low values of the
threshold, there is a range of $\varepsilon$ for which the amount of injury
inflicted to the system is actually higher than the maximum at constant base
damage. This implies that the bad team has the option of causing a higher
injury on a high security system by just engineering the virus with an
internal clock, which would progressively increase the damage inflicted to
infected devices. On the contrary, at high values of $\theta$ (i.e., when
security is not so demanding), viruses whose base damage is given by Eq. (1)
cause an overall injury which is comparable with the maximal value at constant
damage strength.

Finally, we briefly describe the scenario in which the virus is engineered to
increase its base damage not due to an internal clock, but due to a mutation
that occurs all the times the virus is transmitted from an infected device to
a susceptible one. This is tantamount to say that, when a device is infected
at time $t$, the damage caused depends on the story of the specific strain of
the infecting virus. In other words, Eq. (1) is substituted by

$d(t_{i})=\frac{d_{0}e^{\varepsilon t_{i}}}{1+d_{0}(e^{\varepsilon
t_{i}}-1)}\;\;\hbox{with}\;\;d_{0}=0.1,$

where the discrete variable $t_{i}$ is now the number of times that the
infective strain propagated before reaching the actual device. Therefore, the
damage caused is now explicitly a function of the specific history of the
strain propagation. Moreover, for the $I_{u}\to I_{a}$ local transition, the
parameter $\mu$ is determined by looking at the base damage received by the
device at the moment of its contagion.

Figure 5: Mutating virus. $D/N$ vs. $\varepsilon$ (see text for definition),
for different choices of networks (Erdős-Rényi in the left column and
Barabási-Albert in the right column) of mean degree 6, 10 and 14 (see legend
at the bottom of the figure for color code) and different values of the
threshold $\theta$: a,b) $\theta=0.2$, c,d) $\theta=0.6$. For comparison, in
all panels we report also the curves of Fig. 4 (with dashed lines) and the
global maximum of the curves of Fig. 3 (with horizontal dotted lines).

The results are reported in Fig. 5, where $D/N$ is plotted vs. $\varepsilon$
for the case of a mutating virus. One can see that the emerging qualitative
scenario is similar to that of time varying viruses, in that one has an
initial growth and an asymptotic behavior with a local maximum in between.
Moreover, we again observe that for low thresholds a range of $\varepsilon$
values exists for which the injury to the system is higher than the case of a
constant base damage whereas at high thresholds a virus with fixed base damage
is more harmful. The comparison with Fig. 4, however, suggests that a virus
with increasing base damage over time always produces a larger injury to the
system than mutating virus.

In conclusion, we have introduced a novel compartmental model able to describe
the spreading of a malware (and of the awareness of its incidence) on a given
network of devices. The novelty of our approach consists in having considered
vectorial compartments made of two components, the first being descriptive of
the state of the device with respect to the virus propagation, and the second
accounting for the awareness of the device’s user about the presence of the
cyber threat. The model allows to follow the evolution of a cyber-epidemic
from the moment at which a malware is seeded in a network of devices, until
when a given user gets aware of the incurred damage and starts a wave of
awareness which eventually leads to the development of a proper antivirus
software. We then illustrate the overall injury that a malware is able to
produce in Erdős-Rényi and scale-free architectures for both the case in which
the virus is causing a fixed damage on each device and the case in which,
instead, the virus is engineered to mutate while replicating from device to
device. As our world is more and more interconnected, a cyber-epidemic is a
dangerous threat, where a cyber virus would undermine the global Internet
system with catastrophic consequences. Our results are the attempt to
describe, in a customized way, the evolution of such pandemic, and our
conclusion may give hints on how to properly engineer a network of devices to
minimize its vulnerability.

This work has been partially suppported by projects PGC2018-101625-B-I00
(Spanish Ministry, AEI/FEDER, UE) and M1993 (URJC Grant). Authors acknowledge
the usage of the resources, technical expertise and assistance provided by the
supercomputing facility CRESCO of ENEA in Portici (Italy).

## References

  * [1] Ronald Ross.  An application of the theory of probabilities to the study of a priori pathometry.–Part I.  Proceedings of the Royal Society of London. Series A, Containing papers of a mathematical and physical character, 92(638):204–230, 1916.
  * [2] Ronald Ross and Hilda P. Hudson.  An application of the theory of probabilities to the study of a priori pathometry.–Part II.  Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character, 93(650):212–225, 1917.
  * [3] Ronald Ross and Hilda P. Hudson.  An application of the theory of probabilities to the study of a priori pathometry.–Part III.  Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character, 93(650):225–240, 1917.
  * [4] William Ogilvy Kermack, A. G. McKendrick, and Gilbert Thomas Walker.  A contribution to the mathematical theory of epidemics.  Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character, 115(772):700–721, 1927.
  * [5] Stephen Wolfram.  Cellular automata and complexity : collected papers.  Reading, Mass.: Addison-Wesley, 1994.
  * [6] Peter Grassberger.  On the critical behavior of the general epidemic process and dynamical percolation.  Mathematical Biosciences, 63(2):157–172, 1983.
  * [7] Romualdo Pastor-Satorras and Alessandro Vespignani.  Epidemic spreading in scale-free networks.  Physical Review Letters, 86(14):3200, 2001.
  * [8] Stefano Boccaletti, William Ditto, Gabriel Mindlin, and Abdon Atangana.  Modeling and forecasting of epidemic spreading: The case of covid-19 and beyond.  Chaos, Solitons, and Fractals, 135:109794, 2020.
  * [9] Stefano Boccaletti, Gabriel Mindlin, William Ditto, and Abdon Atangana.  Closing editorial: Forecasting of epidemic spreading: lessons learned from the current covid-19 pandemic.  Chaos, Solitons, and Fractals, 139:110278, 2020.
  * [10] Roy M Anderson and Robert M May.  Infectious diseases of humans: dynamics and control.  Oxford university press, 1992.
  * [11] James D Murray.  Mathematical biology. Second corrected edition.  Springer-Verlag, Berlin, 1993.
  * [12] Head Written by Jeremy Kaye.  Protecting critical infrastructure from a cyber pandemic.
  * [13] Winfried Gleissner.  A mathematical theory for the spread of computer viruses.  Computers & Security, 8(1):35–41, 1989.
  * [14] Angel Martín del Rey.  Mathematical modeling of the propagation of malware: a review.  Security and Communication Networks, 8(15):2561–2579, 2015.
  * [15] Soodeh Hosseini and Mohammad Abdollahi Azgomi.  The dynamics of an seirs-qv malware propagation model in heterogeneous networks.  Physica A: Statistical Mechanics and its Applications, 512:803–817, 2018.
  * [16] JD Hernández Guillén and A Martín del Rey.  Modeling malware propagation using a carrier compartment.  Communications in Nonlinear Science and Numerical Simulation, 56:217–226, 2018.
  * [17] Soodeh Hosseini and Mohammad Abdollahi Azgomi.  The dynamics of an seirs-qv malware propagation model in heterogeneous networks.  Physica A: Statistical Mechanics and its Applications, 512:803–817, 2018.
  * [18] Wanping Liu and Shouming Zhong.  A novel dynamic model for web malware spreading over scale-free networks.  Physica A: Statistical Mechanics and Its Applications, 505:848–863, 2018.
  * [19] Asma Al Kindi, Dawood Al Abri, Ahmed Al Maashri, and Fahad Bait-Shiginah.  Analysis of malware propagation behavior in social internet of things.  International Journal of Communication Systems, 32(15):e4102, 2019\.
  * [20] Asma Al Kindi, Dawood Al Abri, Ahmed Al Maashri, and Fahad Bait-Shiginah.  Analysis of malware propagation behavior in social internet of things.  International Journal of Communication Systems, 32(15):e4102, 2019\.
  * [21] Clara Granell, Sergio Gómez, and Alex Arenas.  Dynamical interplay between awareness and epidemic spreading in multiplex networks.  Physical Review Letters, 111(12):128701, 2013.
  * [22] Our model has been simulated using Python, and in particular the existing NetworkX library, with transition functions from one state to the other ruled by the respective probabilities. Unless otherwise specified, the reported numerical results refer to an ensemble average over 500 different network realizations of $n=1000$ nodes each, and a further ensemble average over 20 initial conditions per network (i.e., 20 different random positions where the initial seed of infected-unaware is placed).
  * [23] Albert-László Barabási and Réka Albert.  Emergence of scaling in random networks.  Science, 286(5439):509–512, 1999.
  * [24] Réka Albert and Albert-László Barabási.  Statistical mechanics of complex networks.  Reviews of Modern Physics, 74(1):47, 2002.
  * [25] Paul Erdös and Alfred Rényi.  On random graphs I.  Publicationes Mathematicae Debrecen, 6:290–297, 1959.
  * [26] Réka Albert, Hawoong Jeong, and Albert-László Barabási.  Error and attack tolerance of complex networks.  Nature, 406(6794):378–382, 2000.


1Center for Intelligent & Interactive Robotics Research, Korea Institute of
Science and Technology,
Seoul, 02792, Korea. Email: (614503@kist.re.kr<EMAIL_ADDRESS>∗
Corresponding author

††thanks:  This work was supported by Korea Institute of Science and
Technology (KIST), under Grant 2E32302.

# Dimensionally Homogeneous Jacobian using Extended Selection Matrix for
Performance Evaluation and Optimization of Parallel Manipulators

Hassen Nigatu1 and Doik Kim1∗

###### Abstract

This paper proposes a new methodology for deriving a point-based dimensionally
homogeneous Jacobian, intended for performance evaluation and optimization of
parallel manipulators with mixed degrees of freedom. Optimal manipulator often
rely on performance indices obtained from the Jacobian matrix. However, when
manipulators exhibit mixed translational and rotational freedoms, the
conventional Jacobian’s inconsistency of units lead to unbalanced optimal
result. Addressing this issue, a point-based dimensionally homogeneous
Jacobian has appeared as a prominent solution. However, existing point-based
approaches for formulating dimensionally homogeneous Jacobian are applicable
to a limited variety of parallel manipulators. Moreover, they are complicated
and less intuitive. This paper introduces an extended selection matrix that
combines component velocities from different points to describe the entire
motion of moving plate. This proposed approach enables us to formulate an
intuitive point-based, dimensionally homogeneous Jacobian, which can be
applied to a wide variety of constrained parallel manipulators. To prove the
validity of proposed method, a numerical example is provided utilizing a four-
degree-of-freedom parallel manipulator.

###### keywords:

Dimensionally homogeneous Jacobian, Selection matrix, Inverse Jacobian,
Parallel manipulator, Performance indices

##  1 Introduction

Performance evaluation and obtaining optimized architectural parameters are a
vital step in the design of parallel manipulators (PMs), as these
significantly influences the effectiveness and accuracy of a robot’s
movements. The challenge lies in performing the tasks when the manipulators
degree of freedom (DoFs) is a combination of rotational and translational
types. This is primarily due to the inconsistency in the units or dimensions
of the Jacobian, a factor that significantly affects the performance measuring
indices of parallel manipulators [1, 2].

Several approaches have been suggested to address this problem [3, 4], and
among them, Jacobian-based methods have been widely used. This popularity can
be attributed to their capability to effectively translate the inherent
mapping process from joint velocities to end-effector velocities, providing a
good intuitive framework [1, 5, 6]. There are also variety of Jacobian-based
approaches of homogenizing the units of the Jacobian matrix [4, 7]. Among
these approaches, point-based approach is more intuitive [6, 8]. Despite its
advantage, the firsly introduced point-based dimensionally homogeneous
Jacobian (DHJ) formulations comprise dependent motions in their entries,
resulting in a condition number with unclear physical meaning and potential
erroneous results [6].

In response to this problem, Pond et al. [6] proposed a method to eliminate
the undesired dependent motions from the system. However, the method is quite
complicated to comprehend and involved a tedious derivative procedures which
leads to higher computation cost [8]. To overcome this issue, the selection
matrix with the shifting property and conventional Jacobian is used to
formulate a point-based dimensionally homogeneous Jacobian matrix [8].
However, this previous paper by the authors focused on a specific scenario
where each component’s velocity encompassed the desired motion of the moving
plate, such as 1T2R PMs with $TzRxRy$ type of motion.
Considering the aforementioned limitations, this paper formulates an $f\times
f$, with $f$ representing the DoF of the mechanism, point-based DHJ matrix
using an extended selection matrix. This Jacobian matrix maps the platform’s
nominal linear velocity to the joint rate. Here, nominal linear velocity
refers to the velocity obtained by combining component velocities from
different points, which can represent the entire motion of the moving plate.
This approach integrates the extended selection matrix, the linear velocity of
points on the moving plate and the manipulator’s conventional Jacobian,
resulting in a square DHJ. The dimensionally homogeneity of the resulting
Jacobian is analytically proven. To validate the correctness of the proposed
method, a numerical comparison is carried out using a four-degree-of-freedom
parallel manipulator as an example. First, the distribution of the condition
number is evaluated across the manipulator’s rotational workspace,
highlighting the disparity in the condition number values of the conventional
and dimensionally homogeneous Jacobian. Then, the unit of geometric parameters
are changed from millimeters to meters, and the condition number is reassessed
to determine if it is invariant under unit changes.

##  2 Formulation of the Dimensionally Homogeneous Jacobian

The derivation method of DHJ involves the following steps. First, the screw-
based constraint embedded inverse Jacobian of the manipulator is formulated
and inverted to get the constraint compatible forward relation. Then, points
that might adequately represent the motion of the moving plate are chosen and
related to the Cartesian velocity. Next, the extended selection matrix is
derived and applied to the points’ linear velocity. This will combine
components from different points to effectively describe the moving plate’s
motion, while also eliminating unwanted or dependent components from the
equation. The resulting velocity is termed as the nominal linear velocity.
Finally, the nominal linear velocity of the moving plate and the forward
velocity equation are related with an $f\times f$ dimensionally homogeneous
Jacobian matrix.

###  2.1 Constraint-Embedded Velocity Relation

The screw-based Jacobian of the manipulator can be analytically obtained using
the method introduced in [5, 9, 10, 11]. Given the task velocity,
$\boldsymbol{\dot{\mathscr{x}}}$, of the moving moving plate, the general
inverse velocity equation of the parallel manipulator has the following form.

$\begin{bmatrix}\dot{\boldsymbol{q}}\\\
\boldsymbol{0}\end{bmatrix}=\begin{bmatrix}\boldsymbol{G}_{a}^{T}\\\
\boldsymbol{G}_{c}^{T}\end{bmatrix}\boldsymbol{\dot{\mathscr{x}}}=\begin{bmatrix}\boldsymbol{G}_{av}^{T}&\boldsymbol{G}_{aw}^{T}\\\
\boldsymbol{G}_{cv}^{T}&\boldsymbol{G}_{cw}^{T}\end{bmatrix}\begin{bmatrix}\boldsymbol{v}\\\
\boldsymbol{\omega}\end{bmatrix}$ (1)

The units of the entries in $\boldsymbol{G}$ in Eq. (1), are dependent on the
type of actuators employed in the manipulator. This paper focus exclusively on
scenarios where the manipulator employs only linear or rotational actuators,
not considering situations involving a combination of these actuator types.

Inverting Eq. (1) yields a constraint compatible forward velocity relation as

$\boldsymbol{\dot{\mathscr{x}}}=\boldsymbol{J}\dot{\boldsymbol{q}}=\begin{bmatrix}\boldsymbol{J}_{a}&\boldsymbol{J}_{c}\end{bmatrix}\begin{bmatrix}\dot{\boldsymbol{q}}_{a}\\\
\boldsymbol{0}\end{bmatrix}$ (2)

In Eq. (2), $\boldsymbol{J}\in{\mathbb{R}}^{6\times 6}$ is the inverse of
$\boldsymbol{G}^{T}$ and its sub-matrix $\boldsymbol{J}_{c}$ is related to the
constraint. Thus, we can explicitly describe the relation of
$\boldsymbol{\dot{\mathscr{x}}}$ and $\dot{\boldsymbol{q}}_{a}$ as

$\boldsymbol{\dot{\mathscr{x}}}=\boldsymbol{J}_{a}\dot{\boldsymbol{q}}_{a},~{}\text{where}~{}\boldsymbol{J}_{a}=\begin{bmatrix}\boldsymbol{J}_{a1}\\\
\boldsymbol{J}_{a2}\end{bmatrix}$ (3)

The Cartesian velocity,
$\boldsymbol{\dot{\mathscr{x}}}\in{\mathbb{R}}^{6\times 1}$, in Eq. (3) is
constraint compatible.When the manipulator employs linear actuators,
$\boldsymbol{J}_{a1}$ is dimensionless, while $\boldsymbol{J}_{a2}$ has a unit
of $\frac{1}{\text{length}}$. Conversely, if the manipulator utilizes
rotational actuators, $\boldsymbol{J}_{a1}$ has a unit of length and
$\boldsymbol{J}_{a2}$ dimensionless. Considering these distinctions, the
point’s linear velocity and selection matrix are established to ensure
consistency or removal of units in the Jacobian.

###  2.2 Linear Velocity of Points

According to the well known shifting property [12] in the rigid body
kinematics, any points velocity on the moving plate can be related to the
Cartesian velocity of the moving plate as

$\boldsymbol{v}_{i}=\boldsymbol{v}+\boldsymbol{\omega}\times\boldsymbol{a}_{i}$
(4)

where $\boldsymbol{v}$ and $\boldsymbol{\omega}$ denotes the linear and
angular velocity of the moving plate, while $\boldsymbol{a}_{i}$ corresponds
to a constant vector extending from the origin of the Cartesian reference
frame to the $i^{th}$ point on the moving plate. Expanding Eq. (4) reveals the
motion of the moving plate that each component of $\boldsymbol{v}_{i}$
encompasses.

$\displaystyle{v}_{ix}$
$\displaystyle={v}_{x}+{\omega}_{y}{a}_{iz}-{\omega}_{z}{a}_{iy}$ (5)
$\displaystyle{v}_{iy}$
$\displaystyle={v}_{y}-{\omega}_{x}{a}_{iz}+{\omega}_{z}{a}_{ix}$
$\displaystyle{v}_{iz}$
$\displaystyle={v}_{z}+{\omega}_{x}{a}_{iy}-{\omega}_{y}{a}_{ix}$

By distributing these points on the moving plate in a noncollinear manner, it
is possible to satisfy the minimum requirement of points needed to fully
represent the motion of the moving plate. Theoretically, the translations of
three noncollinear points on the moving plate are sufficient to uniquely
identify the motion of the body in terms of translation and rotation, but more
points may be required depending on the DoF of the mechanism.

Hence, Eq. (4) can be generalized as

$\displaystyle\boldsymbol{v}_{p}=\begin{bmatrix}\boldsymbol{v}_{1}\\\
\vdots\\\ \boldsymbol{v}_{i}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}\boldsymbol{I}&-[\boldsymbol{a}_{1}]_{\times}\\\
\vdots&\vdots\\\
\boldsymbol{I}&-[\boldsymbol{a}_{i}]_{\times}\end{bmatrix}\begin{bmatrix}\boldsymbol{v}\\\
\boldsymbol{\omega}\end{bmatrix}$ (6)
$\displaystyle=\boldsymbol{V}_{p}\boldsymbol{\dot{\mathscr{x}}}$

where $\boldsymbol{V}_{p}\in{\mathbb{R}}^{3f\times 6}$ maps the moving plate
cartesian velocity to the points velocity on the moving plate. Vector
$\boldsymbol{v}_{i}$ in Eq. (6) has three components and hence from
$\boldsymbol{v}_{p}\in{\mathbb{R}}^{3f\times 1}$, we need to determine the
components that can appropriately describe the motion of the moving plate via
a selection matrix [8] as follows

$\displaystyle\boldsymbol{S}\boldsymbol{v}_{p}$
$\displaystyle=\boldsymbol{S}\boldsymbol{V}_{p}\boldsymbol{\dot{\mathscr{x}}},~{}\text{where}~{}\boldsymbol{S}\in{\mathbb{R}}^{f\times
3f}~{}\text{is a selection}$ (7) $\displaystyle\text{matrix
that}~{}\text{extracts the components from}~{}\boldsymbol{v}_{p}.$
$\displaystyle\boldsymbol{v}_{ps}$
$\displaystyle=\boldsymbol{V}_{ps}\boldsymbol{\dot{\mathscr{x}}}~{}\text{where}~{}\boldsymbol{V}_{ps}\in{\mathbb{R}}^{f\times
6}$

However, deriving the selection matrix $\boldsymbol{S}$ is not always
straightforward. This is because only manipulators whose moving plates exhibit
$T_{x}R_{y}R_{z}$, $T_{y}R_{x}R_{z}$ and $T_{z}R_{x}R_{y}$ types of motion can
be uniquely represented with the component velocities shown in Eq. (5). For a
comprehensive understanding of the establishment of selection matrices for
these groups of PMs, readers are encouraged to refer to [8]. PMs falling
outside of these categories will need to utilize a combination of components
from different points, an approach that is covered in this paper.

###  2.3 Dimensionally Homogeneous Jacobian

In this paper, we derive the dimensionally homogeneous Jacobian by
representing the motion of the moving plate using linear velocity, ensuring
uniform units across its entries. However, it is important to note that the
linear velocities used here are not merely the component velocities of
individual points on the moving plate. Instead, they are a combination of
components from various points. This approach is used to encompass all desired
motion of the moving plate into a representative velocity equation, which we
call it the nominal velocity.

To derive the dimensionally homogeneous Jacobian, relations, Eq. (3) and Eq.
(7) are combined as follows

$\displaystyle\boldsymbol{v}_{ps}$
$\displaystyle=\boldsymbol{V}_{ps}\boldsymbol{\dot{\mathscr{x}}}$ (8)
$\displaystyle=\boldsymbol{V}_{ps}\boldsymbol{J}_{a}\dot{\boldsymbol{q}}_{a}$
$\displaystyle=\boldsymbol{J}_{dh}\dot{\boldsymbol{q}}_{a}$

In Eq. (8), $\boldsymbol{J}_{dh}\in{\mathbb{R}}^{f\times f}$ is a Jacobian
that relates the nominal linear velocity ($\boldsymbol{v}_{ps}$) of the moving
plate to the actuated joint rate
($\dot{\boldsymbol{q}}_{a}\in{\mathbb{R}}^{f\times 1}$). To demonstrate the
consistency of units in its entries, we considered the following two generic
cases.

Case 1: PMs with linear actuators. In this case, $\dot{\boldsymbol{q}}_{a}$
has unit of $\frac{\text{length}}{\text{time}}$ while $\boldsymbol{S}$ is
dimensionless. Referring to Eq. (6), we can observe that the first term is
dimensionless, while the second term has a unit of length. Furthermore, in Eq.
(2), we know Block matrix $\boldsymbol{J}_{a1}$ is dimensionless and
$\boldsymbol{J}_{a2}$ has a unit of $\frac{\text{1}}{\text{length}}$. As a
result, we conclude that the Jacobian for this particular group of
manipulators is dimensionless.

Case 2: PMs with rotational actuators. For PMs with rotational actuators
$\dot{\boldsymbol{q}}_{a}$ has unit of $\frac{\text{angle}}{\text{time}}$
while the unit of $\boldsymbol{V}_{p}$ is unchanged. Furthermore, the matrix
$\boldsymbol{J}_{a1}$ has a unit of length and $\boldsymbol{J}_{a2}$ is
dimensionless for this group of PMs. Consequently, the resulting Jacobian
$\boldsymbol{J}_{dh}$ to has a unit of length which is consistent.

Because entries of $\boldsymbol{J}_{dh}$ are either dimensionless or
dimensionally homogeneous, its condition number or singular values have
physical significance and can be used to measure the dexterity of the
manipulator.

The next section demonstrates how to derive it by considering a relevant
example: a four DoF (degrees of freedom) $T_{y}T_{z}R_{x}R_{y}$ type Parallel
Manipulator (PM).

##  3 Example

The mechanism depicted in Fig. 1 is a $T_{y}T_{z}R_{x}R_{y}$ type 4 DoF PM
[13] with a PUS joint order in the first and third limbs, and a PRS type joint
sequence in the second and fourth limbs. The P joint is parallel to the
$z$-axis, while the R joint in the PRS limb is parallel to the $x$-axis, and
the U joint in the PUS limb has axes parallel to the $x$ and $y$ axes,
respectively. The mechanism is capable of rotating about the $x$ and $y$ axes,
as well as translating along the $y$ and $z$ directions. However, due to the
presence of revolute joints in the second and fourth limbs, the mechanism is
constrained in terms of $x$-axis translation and $z$-axis rotations, making it
a zero-torsion type PM mechanism. The DoF of the mechanism can also be
determined by employing Tsai’s DoF formula, which is expressed as follows:

$\displaystyle F$ $\displaystyle=\lambda(n-j-1)+\sum_{i=1}^{j}f_{i}$ (9)
$\displaystyle 4$ $\displaystyle=6(10-12-1)+22$

Point $A_{i}$ at the based is location of limbs while $B_{i}$ is the center of
spherical joints. Point $C_{i}$ is the center universal joints for the first
and third limbs while it is the center of revolute joints for limbs 2 and 4.
The position vector $\boldsymbol{a}_{i}$ is extended from origin moving frame
to the $i^{th}$ spherical joint while $\boldsymbol{b}_{i}$ is extended from
the fixed frame to the point $A_{i}$. The direction vector
$\boldsymbol{s}_{ji\parallel}$ is associated to each joint axis.

Figure 1: Four DoF PM

To appropriately represent the motion of the moving plate, we need four points
and for convenience these points are chosen to be the center of the spherical
joints. Expanding Eq. (4) to the four points located at the center of
spherical joints at the moving plate, we get a $12\times 6$ matrix that
relates the $i^{th}$ points linear velocity to the moving plate center
velocity ($\boldsymbol{\dot{\mathscr{x}}}$) as

$\begin{bmatrix}v_{1x}\\\ v_{1y}\\\ v_{1z}\\\ \vdots\\\ v_{4x}\\\ v_{4y}\\\
v_{4z}\end{bmatrix}=\begin{bmatrix}1&0&0&0&a_{1z}&-a_{1y}\\\
0&1&0&-a_{1z}&0&a_{1x}\\\ 0&0&1&a_{1y}&-a_{1x}&0\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\\ 1&0&0&0&a_{4z}&-a_{4y}\\\
0&1&0&-a_{4z}&0&a_{4x}\\\
0&0&1&a_{4y}&-a_{4x}&0\end{bmatrix}\begin{bmatrix}v_{x}\\\ v_{y}\\\ v_{z}\\\
\omega_{x}\\\ \omega_{y}\\\ \omega_{z}\end{bmatrix}$ (10)

The points’ linear velocity in Eq. (10) includes 12 components, three for each
point. Therefore, we have many dependent motions, yet we only require four
components. Referring to Eq. (5), none of the components encompass the desired
motion of the moving plate. Hence, we need to formulate a selection matrix
that can combine components from different points and obtain a nominal
velocity that describes the motion of the moving plate. As the independent
motion of the moving plate for this manipulator are $v_{y},v_{z},\omega_{x}$
and $\omega_{y}$, combining $v_{iy}$ and $v_{iz}$ components can sufficiently
describe the manipulator’s motion. However, a combination of component is not
unique and one can freely choose one of the following pairs.

Limb 1: $\displaystyle~{}(v_{1y},v_{2z}),(v_{1y},v_{3z}),(v_{1y},v_{4z})$ (11)
Limb 2: $\displaystyle~{}(v_{2y},v_{1z}),(v_{2y},v_{3z}),(v_{2y},v_{4z})$ Limb
3: $\displaystyle~{}(v_{3y},v_{1z}),(v_{3y},v_{2z}),(v_{3y},v_{4z})$ Limb 4:
$\displaystyle~{}(v_{4y},v_{1z}),(v_{4y},v_{2z}),(v_{4y},v_{3z})$

For this particularly case, we selected the following combination from Eq.
(11).

Limb 1: $\displaystyle~{}(v_{1y},v_{2z})$ (12) Limb 2:
$\displaystyle~{}(v_{2y},v_{3z})$ Limb 3: $\displaystyle~{}(v_{3y},v_{4z})$
Limb 4: $\displaystyle~{}(v_{4y},v_{1z})$

By utilizing Eq. (12), we can establish the extended selection matrix as

$\displaystyle\boldsymbol{S}$
$\displaystyle=\left[\begin{matrix}0&-\dfrac{a_{2x}}{a_{1x}-a_{2x}}&1&0&\dfrac{a_{1x}}{a_{1x}-a_{2x}}&0\\\
0&0&0&0&-\dfrac{a_{3x}}{a_{2x}-a_{3x}}&1\\\ 0&0&0&0&0&0\\\
0&-\dfrac{a_{4x}}{a_{1x}-a_{4x}}&1&0&0&0\\\ \end{matrix}\right.\ $ (13)
$\displaystyle\qquad\ \left.\begin{matrix}0&0&0&0&0&0\\\
0&\dfrac{a_{2x}}{a_{2x}-a_{3x}}&0&0&0&0\\\
0&-\dfrac{a_{4x}}{a_{3x}-a_{4x}}&1&0&\dfrac{a_{3x}}{a_{3x}-a_{4x}}&0\\\
0&0&0&0&\dfrac{a_{1x}}{a_{1x}-a_{4x}}&0\end{matrix}\right]$

It is quite important to note that $\boldsymbol{S}$ in Eq. (13) is
dimensionless and same as to that the usual selection matrix in terms of
units. However, the usual selection matrices [8] have entries of 1s and 0s
unlike the extended selection matrix derived in the paper. Then, multiplying
Eq. (13) with Eq. (10), matrix$\boldsymbol{V}_{ps}\in{\mathbb{R}}^{4\times 4}$
relates the nominal velocity ($\boldsymbol{v}_{ps}$) and independent Cartesian
velocity of the moving plate as in Eq. (14).

$\begin{aligned} \begin{bmatrix}\underline{v_{1}}\\\ \underline{v_{2}}\\\
\underline{v_{3}}\\\
\underline{v_{4}}\end{bmatrix}&=\left[\begin{matrix}\dfrac{a_{1x}}{a_{1x}-a_{2x}}-\dfrac{a_{2x}}{a_{1x}-a_{2x}}&1\\\
\dfrac{a_{2x}}{a_{2x}-a_{3x}}-\dfrac{a_{3x}}{a_{2x}-a_{3x}}&1\\\
\dfrac{a_{3x}}{a_{3x}-a_{4x}}-\dfrac{a_{4x}}{a_{3x}-a_{4x}}&1\\\
\dfrac{a_{1x}}{a_{1x}-a_{4x}}-\dfrac{a_{4x}}{a_{1x}-a_{4x}}&1\\\
\end{matrix}\right.\ \\\ &\qquad\
\left.\begin{matrix}\dfrac{a_{1y}a_{1x}-a_{2x}a_{1z}+a_{1x}a_{3z}}{a_{1x}-a_{2x}}&-a_{1x}\\\
\dfrac{a_{3y}a_{2x}+a_{2x}a_{3z}-a_{3x}a_{3z}}{a_{2x}-a_{3x}}&-a_{2x}\\\
\dfrac{a_{3y}a_{3x}+a_{3x}a_{4z}-a_{4x}a_{3z}}{a_{3x}-a_{4x}}&-a_{3x}\\\
\dfrac{a_{1y}a_{1x}+a_{1x}a_{4z}-a_{4x}a_{1z}}{a_{1x}-a_{4x}}&-a_{1x}\end{matrix}\right]\end{aligned}\begin{bmatrix}vy\\\
vz\\\ wx\\\ wy\end{bmatrix}$ (14)

where, $\begin{bmatrix}\underline{v_{1}}\\\ \underline{v_{2}}\\\
\underline{v_{3}}\\\
\underline{v_{4}}\end{bmatrix}=\begin{bmatrix}\dfrac{(v_{2y}+v_{1z})a_{1x}-(v_{1y}+v_{1z})a_{2x}}{a_{1x}-a_{2x}}\\\
\dfrac{(v_{3y}+v_{2z})a_{2x}-(v_{2y}+v_{2z})a_{3x}}{a_{2x}-a_{3x}}\\\
\dfrac{(v_{4y}+v_{3z})a_{3x}-(v_{3y}+v_{3z})a_{4x}}{a_{3x}-a_{4x}}\\\
\dfrac{(v_{4y}+v_{1z})a_{1x}-(v_{1y}+v_{1z})a_{4x}}{a_{1x}-a_{4x}}\\\
\end{bmatrix}$

The inverse Jacobian of the manipulator is obtained through the analytic screw
theory method and is given as shown in Eq. (15). The first four row of
$\boldsymbol{G}^{T}$ represents the motion Jacobian while the last two depicts
the structural constraints. Hence,
$\boldsymbol{G}_{c}^{T}\boldsymbol{\dot{\mathscr{x}}}=\boldsymbol{0}$ is
always satisfied.

$\dot{\boldsymbol{q}}=\begin{bmatrix}\boldsymbol{G}_{a}^{T}\\\
\boldsymbol{G}_{c}^{T}\end{bmatrix}\boldsymbol{\dot{\mathscr{x}}}=\begin{bmatrix}\dfrac{\boldsymbol{n}_{1}^{T}}{\boldsymbol{n}_{1}^{T}\boldsymbol{s}_{11\parallel}}&\dfrac{(\boldsymbol{n}_{1}\times\boldsymbol{a}_{1})^{T}}{\boldsymbol{n}_{1}^{T}\boldsymbol{s}_{11\parallel}}\\\
\dfrac{\boldsymbol{l}_{2}^{T}}{\boldsymbol{l}_{2}^{T}\boldsymbol{s}_{12\parallel}}&\dfrac{(\boldsymbol{l}_{2}\times\boldsymbol{a}_{2})^{T}}{\boldsymbol{l}_{2}^{T}\boldsymbol{s}_{12\parallel}}\\\
\dfrac{\boldsymbol{n}_{3}^{T}}{\boldsymbol{n}_{3}^{T}\boldsymbol{s}_{13\parallel}}&\dfrac{(\boldsymbol{n}_{3}\times\boldsymbol{a}_{3})^{T}}{\boldsymbol{n}_{3}^{T}\boldsymbol{s}_{13\parallel}}\\\
\dfrac{\boldsymbol{l}_{4}^{T}}{\boldsymbol{l}_{4}^{T}\boldsymbol{s}_{14\parallel}}&\dfrac{(\boldsymbol{l}_{4}\times\boldsymbol{a}_{4})^{T}}{\boldsymbol{l}_{4}^{T}\boldsymbol{s}_{14\parallel}}\\\
\boldsymbol{s}_{22\parallel}^{T}&(\boldsymbol{s}_{22\parallel}\times\boldsymbol{a}_{2})^{T}\\\
\boldsymbol{s}_{24\parallel}^{T}&(\boldsymbol{s}_{24\parallel}\times\boldsymbol{a}_{4})^{T}\\\
\end{bmatrix}\begin{bmatrix}\boldsymbol{v}\\\
\boldsymbol{\omega}\end{bmatrix}$ (15)

where $\boldsymbol{G}_{a}^{T}\in{\mathbb{R}}^{4\times 6}$ and
$\boldsymbol{G}_{c}^{T}\in{\mathbb{R}}^{2\times 6}$ and
$\boldsymbol{n}_{i}=\boldsymbol{s}_{3i\parallel}\times\boldsymbol{s}_{2i\parallel}$.
$\boldsymbol{l}_{i}$ is a vector extending from $C_{i}$ to $B_{i}$.

The first term of $\boldsymbol{G}^{T}$ is dimensionless while the second term
has a unit of length. Hence, units of the inverse Jacobian of this manipulator
is inconsistent and must be changed to dimensionless or consistent unit.

The forward Jacobian $\boldsymbol{J}_{a}\in{\mathbb{R}}^{6\times 3}$ is
analytically obtained by inverting $\boldsymbol{G}^{-T}$ as in Eq. (16).

$\boldsymbol{J}_{a}=\\!\begin{aligned}
&\left[\begin{matrix}\boldsymbol{G}_{av}^{-T}\boldsymbol{G}_{aw}^{T}(\boldsymbol{G}_{cw}^{T}-\boldsymbol{G}_{cv}^{T}\boldsymbol{G}_{av}^{-T}\times\\\
-(\boldsymbol{G}_{cw}^{T}-\boldsymbol{G}_{cv}^{T}\boldsymbol{G}_{av}^{-T}\boldsymbol{G}_{aw})^{-1}\times\end{matrix}\right.\\\
&\qquad\qquad\left.\begin{matrix}{}\boldsymbol{G}_{av}^{-T}+\boldsymbol{G}_{aw}^{T})^{-1}\boldsymbol{G}_{cv}^{T}\boldsymbol{G}_{av}^{-T}\\\
{}\boldsymbol{G}_{cv}^{T}\boldsymbol{G}_{av}^{-T}\end{matrix}\right]\end{aligned}$
(16)

By substituting Eq. (13) and Eq. (10) into Eq. (7), and subsequently replacing
$\boldsymbol{\dot{\mathscr{x}}}$ with $\boldsymbol{J}_{a}\dot{\boldsymbol{q}}$
in Eq. (7), we derive the $4\times 4$ dimensionless Jacobian as discussed in
case 1.

###  3.1 Numerical Evaluation

In order to verify the correctness of the derived dimensionally homogeneous
Jacobian, the distribution of the condition number ($k$) for the manipulator
over the entire workspace is evaluated using geometric and motion parameters
outlined in Table 1.

Table 1: Structural and pose parameters of the PM. Parameters | value
---|---
Radius of the moving plate $(r_{a})$ | 200 mm
Radius of the base plate $(r_{b})$ | 450 mm
Fixed length link $(l)$ | 687 mm
$\theta$ | $\pm 50^{o}$
$\psi$ | $\pm 50^{o}$
$y$ | 0 mm
$z$ | 100-200 mm

It is known that in parallel manipulator design, the condition number ($k$) of
the Jacobian matrix can be used as a performance measure to evaluate the
quality of motion, precision, and stability of the manipulator. The best value
of $k$ is 1 which is the minimum possible value and it indicates that all
columns (or rows) of the Jacobian matrix matrix are orthogonal to each other.
This implies that the system of equations is well-conditioned and the solution
will not be overly sensitive to errors in the data or to small changes in the
input. This can be interpreted as that the manipulator is isotropic [2]. As
$k$ increases beyond 1, the system of equations becomes increasingly ill-
conditioned. This means that the solution may be very sensitive to errors in
the data or to small changes in the input and hence the manipulator is
approaching to singularity. Contrary to this, if $k$ values is small enough or
remains closer to 1, it can be interpreted the manipulator is away from
singular configuration.

Accordingly, $k=cond(\boldsymbol{G}^{T})$ is first computed and the result is
shown in Fig. 2 over the rotational workspace. The simulation result indicate
a substantial increase in the condition number, which do not to adequately
reflect the physical properties of the manipulator. Consequently,
$cond(\boldsymbol{J}_{dh})$ was determined, as shown in Fig. 3. In the
rotational workspace, the value of $k$ remained low and near 1. This value of
$k$ can properly indicate whether the manipulator is far or approaching a
singular configuration.

Figure 2: $k=cond(\boldsymbol{G})$ Figure 3: $k=cond(\boldsymbol{J}_{dh})$
Figure 4: $k=cond(\boldsymbol{J}_{dh})$

Additionally, the sensitivity of $k$ of both Jacobians to unit changes is
evaluated over the rotational workspace. The results shows significant
discrepancy in the condition number of the conventional Jacobian when units
are changed from millimeter to meters, as depicted in Fig. 4. Comparing this
result with Fig. 2 can be considered more optimized even-though nothing is
changed but unit. However, the value of $cond(\boldsymbol{J}_{dh})$, when
measured in meters, remained unchanged. The result of
$cond(\boldsymbol{J}_{dh})$ in meter is not provided here because it is the
same as to that of shown in Fig. 3. Hence, $cond(\boldsymbol{J}_{dh})$ is
invariant under the change of units. As previously mentioned, the choice of
component combinations is not unique. Hence, we have the flexibility to choose
various pairs of $v_{iy}$ and $v_{iz}$ from the provided candidates. For
instance, by selecting $(v_{1y},v_{3z}),(v_{2y},v_{4z}),(v_{3y},v_{1z})$, and
$(v_{4y},v_{2z})$, we can derive the following selection matrix.

$\displaystyle\boldsymbol{S}$
$\displaystyle=\left[\begin{matrix}0&-\dfrac{a_{2x}}{a_{3x}-a_{2x}}&1&0&\dfrac{a_{3x}}{a_{3x}-a_{2x}}&0\\\
0&0&0&0&-\dfrac{a_{4x}}{a_{2x}-a_{4x}}&1\\\ 0&0&0&0&0&0\\\
0&-\dfrac{a_{4x}}{a_{2x}-a_{4x}}&1&0&0&0\\\ \end{matrix}\right.\ $
$\displaystyle\qquad\ \left.\begin{matrix}0&0&0&0&0&0\\\
0&\dfrac{a_{2x}}{a_{2x}-a_{4x}}&0&0&0&0\\\
0&-\dfrac{a_{1x}}{a_{3x}-a_{1x}}&1&0&\dfrac{a_{3x}}{a_{3x}-a_{1x}}&0\\\
0&0&0&0&\dfrac{a_{2x}}{a_{2x}-a_{4x}}&0\end{matrix}\right]$

Then, with this selection matrix, we establish the dimensionally homogeneous
Jacobian as shown in Eq. (8) and the condition number distribution over the
workspace is evaluated. The simulation has shown the same result to that of
the dimensionally homogeneous Jacobian obtained using Eq. (13).

This consistent property of the Jacobian matrix is quite important while using
the condition number as a measure of performance or computing the dexterity
for parameter optimization of PMs.

##  4 Conclusion

This paper introduces an extended selection matrix to formulate a point-based,
dimensionally homogeneous Jacobian of various constrained parallel
manipulators. The proposed method allows for the derived Jacobian’s condition
number and singular values to be utilized as a performance index and
optimization with unit independence.

To validate the proposed approach, the condition number ($k$) for both the
conventional Jacobian ($\boldsymbol{G}$) and the dimensionally homogeneous
Jacobian ($\boldsymbol{J}_{dh}$) across the rotational workspace were
compared. Simulation results indicated a large value of $k$ for
$\boldsymbol{G}$ and a remarkably stable value of $k$ for
$\boldsymbol{J}_{dh}$.

Further, we reassessed the distribution of the $k$ value for the two Jacobians
by changing the units from millimeters to meters. The results confirmed that
$k$ of $\boldsymbol{G}$ varied significantly, while $k$ of
$\boldsymbol{J}_{dh}$ remained consistent, irrespective of the unit change.
This phenomenon proves the dimensional homogeneity of the proposed Jacobian,
where both the linear and angular parts exhibit similar value distributions
and are not unit-dependent. As a result, our method allows for the correct
optimization of the manipulators with mixed DoFs. By employing the proposed
approach for different manipulators with mixed DoFs, we can confidently assess
and optimize their performance.

## ACKNOWLEDGMENT

This work was supported by Korea Institute of Science and Technology (KIST),
under Grant 2E32302.

## References

  * [1] C. M. Gosselin, “The optimum design of robotic manipulators using dexterity indices,” _Robotics and Autonomous Systems_ , vol. 9, no. 4, pp. 213–226, 1992.
  * [2] J. P. Merlet, “Jacobian, Manipulability, Condition Number and Accuracy of Parallel Robots,” in _Robotics Research_. Berlin, Heidelberg: Springer Berlin Heidelberg, 2007, vol. 28, pp. 175–184.
  * [3] C. Yang, W. Ye, and Q. Li, “Review of the performance optimization of parallel manipulators,” _Mechanism and Machine Theory_ , vol. 170, no. November 2021, p. 104725, 2022.
  * [4] A. Rosyid, B. El-Khasawneh, and A. Alazzam, “Review article: Performance measures of parallel kinematics manipulators,” _Mechanical Sciences_ , vol. 11, no. 1, pp. 49–73, 2020.
  * [5] “Analytic formulation of reciprocal screws and its application to nonredundant robot manipulators,” _Journal of Mechanical Design, Transactions of the ASME_ , vol. 125, no. 1, pp. 158–164, 2003.
  * [6] G. Pond and J. A. Carretero, “Formulating Jacobian matrices for the dexterity analysis of parallel manipulators,” _Mechanism and Machine Theory_ , vol. 41, no. 12, pp. 1505–1519, dec 2006.
  * [7] H. Nigatu, “Constraint-Embedded Velocity Analysis of Lower-DoF Parallel Manipulators to Remove Parasitic Motion .” Dissertation, UST, 2021.
  * [8] H. Nigatu and D. Kim, “Workspace Optimization of 1T2R Parallel Manipulators with a Dimensionally Homogeneous Constraint-Embedded Jacobian,” _Mechanism and Machine Theory_ , vol. 188, pp. 1–24, 2023.
  * [9] D. Kim, “Kinematic Analysis of Spatial Parallel Manipulators: Analytic Approach with The Restriction Space,” Dissertation, POSTECH, 2002.
  * [10] H. Nigatu, Y. H. Choi, and D. Kim, “On the structural constraint and motion of 3-PRS parallel kinematic machines,” in _Proceedings of the ASME Design Engineering Technical Conference_ , vol. 8A-2021. American Society of Mechanical Engineers (ASME), 2021.
  * [11] H. Nigatu, “Constraint-Embedded Velocity Analysis of Lower-DoF Parallel Manipulators to Remove Parasitic Motion .” Dissertation, UST, 2021. [Online]. Available: https://data.mendeley.com/datasets/pw4dtc2v8r/3
  * [12] H. Lipkin, “Material and local time derivatives of screws with applications to dynamics and stiffness,” _Proceedings of the ASME Design Engineering Technical Conference_ , vol. 2 B, no. 1999, pp. 1339–1348, 2004.
  * [13] Y. Li, E. Zhang, Y. Song, and Z. Feng, “Stiffness modeling and analysis of a novel 4-DOF PKM for manufacturing large components,” _Chinese Journal of Aeronautics_ , vol. 26, no. 6, pp. 1577–1585, dec 2013.


k(G)$ where the equality is by (105) and the inequality is by the definition
of a minimal bounding subset of a bounding group.

  * •

For each bounding group $G\in{\mathcal{B}}({\mathcal{I}})$ it holds that
${\left\lceil\frac{|G\setminus A_{0}|}{k(G)}\right\rceil}\leq
p({\mathcal{I}})-1$.

${\left\lceil\frac{|G\setminus
A_{0}|}{k(G)}\right\rceil}={\left\lceil\frac{|G\setminus
S_{G}|}{k(G)}\right\rceil}\leq p({\mathcal{I}})-1$

The first equality is by (105). The first inequality is because $S_{G}\in
b^{*}(G)$ and in particular $S_{G}\in b(G)$; thus, the inequality follows by
the definition of a bounding subset of $G$.

∎

######  Claim I.2.

For any BPP instance ${\mathcal{I}}=(I,{\mathcal{G}},s,k)$ it holds that
$A_{1}\in{\mathcal{C}}({\mathcal{I}})$ and for each bounding group
$G\in{\mathcal{B}}({\mathcal{I}})$ it holds that
${\left\lceil\frac{|G\setminus A_{1}|}{k(G)}\right\rceil}\leq
p({\mathcal{I}})-1$.

###### Proof.

We prove the claim using loop invariant for the while loop of Step 10. Let
$A(s),A(t)$ be the object $A$ before and after an iteration of the while loop
of Step 10, respectively. Now, assume that
$A(s)\in{\mathcal{C}}({\mathcal{I}})$ and for each bounding group
$G\in{\mathcal{B}}({\mathcal{I}})$ it holds that
${\left\lceil\frac{|G\setminus A(s)|}{k(G)}\right\rceil}\leq
p({\mathcal{I}})-1$. We show the conditions of the claim below for $A(t)$. By
Step 10, let $\ell\in I$ such that $A(t)=A(s)\cup\\{\ell\\}$. Therefore,

  * •

$s(A(t))\leq 1$.

$s(A(t))=s\left(A(s)\right)+s(\ell)\leq(1-\delta)+\delta=1.$

The first equality is by Step 10. The inequality is because $s(A(s))\leq
1-\delta$ by Step 10; in addition, for all items the size is at most $\delta$.

  * •

For all $G\in{\mathcal{G}}$ it holds that $|G\cap A(t)|\leq k(G)$. If
$G\in{\mathcal{G}}\setminus\\{\textsf{group}(\ell)\\}$; then, $|G\cap
A(t)|=|A(s)\cap G|\leq k(G)$. The equality is by Step 10 and the inequality is
by the assumption that $A(s)\in{\mathcal{C}}({\mathcal{I}})$. Otherwise,
$|G\cap A(t)|=\left|A(s)\cap G\right|+|\\{\ell\\}\cap G|=|A(s)\cap G|+1\leq
k(G)$. The first equality is by Step 10 and the inequality is by Step 10.

  * •

For each bounding group $G\in{\mathcal{B}}({\mathcal{I}})$ it holds that
${\left\lceil\frac{|G\setminus A(t)|}{k(G)}\right\rceil}\leq
p({\mathcal{I}})-1$. If
$G\in{\mathcal{G}}\setminus\\{\textsf{group}(\ell)\\}$; then,
${\left\lceil\frac{|G\setminus
A(t)|}{k(G)}\right\rceil}={\left\lceil\frac{|G\setminus
A(s)|}{k(G)}\right\rceil}\leq p({\mathcal{I}})-1$. The equality is by Step 10
and the inequality is by the assumption on $A(s)$. Otherwise,

${\left\lceil\frac{|G\setminus
A(t)|}{k(G)}\right\rceil}={\left\lceil\frac{|G\setminus
A(s)|-|\\{\ell\\}|}{k(G)}\right\rceil}\leq{\left\lceil\frac{|G\setminus
A(s)|}{k(G)}\right\rceil}\leq p({\mathcal{I}})-1.$

The first equality is by Step 10 and the last inequality is by the assumption
on $A(s)$. Therefore, by Claim I.1 and the above loop invariant the claim
follows.

∎

######  Claim I.3.

$A\in{\mathcal{C}}({\mathcal{I}})$ and for each bounding group
$G\in{\mathcal{B}}({\mathcal{I}})$ it holds that
${\left\lceil\frac{|G\setminus A|}{k(G)}\right\rceil}\leq p({\mathcal{I}})-1$.

###### Proof.

We prove the claim using loop invariant for the while loop of Step 10. Let
$A(s),A(t)$ be the object $A$ before and after an iteration of the while loop
of Step 10, respectively. Now, assume that
$A(s)\in{\mathcal{C}}({\mathcal{I}})$ and for each bounding group
$G\in{\mathcal{B}}({\mathcal{I}})$ it holds that
${\left\lceil\frac{|G\setminus A(s)|}{k(G)}\right\rceil}\leq
p({\mathcal{I}})-1$. We show the conditions of the claim below for $A(t)$. By
Step 10, let $\ell,y\in I$ such that
$A(t)=\left(A(s)\setminus\\{\ell\\}\right)\cup\\{y\\}$. Therefore,

  * •

$s(A(t))\leq 1$.

$s(A(t))=s\left(A(s)\right)-s(\ell)+s(y)\leq
s\left(A(s)\right)+s(y)\leq(1-\delta)+\delta=1.$

The first equality is by Step 10. The second inequality is because
$s(A(s))\leq 1-\delta$ by Step 10; in addition, for all items the size is at
most $\delta$.

  * •

For all $G\in{\mathcal{G}}$ it holds that $|G\cap A(t)|\leq k(G)$. If
$G\in{\mathcal{G}}\setminus\\{\textsf{group}(\ell)\\}$; then, $|G\cap
A(t)|=|A(s)\cap G|\leq k(G)$. The equality is by Step 10 and the inequality is
by the assumption that $A(s)\in{\mathcal{C}}({\mathcal{I}})$. Otherwise,
$|G\cap A(t)|=\left|A(s)\cap G\right|-|\\{\ell\\}\cap G|+|\\{y\\}\cap
G|=|A(s)\cap G|\leq k(G)$. The first equality is by Step 10 and the inequality
is by the assumption that $A(s)\in{\mathcal{C}}({\mathcal{I}})$

  * •

For each bounding group $G\in{\mathcal{B}}({\mathcal{I}})$ it holds that
${\left\lceil\frac{|G\setminus A(t)|}{k(G)}\right\rceil}\leq
p({\mathcal{I}})-1$. If
$G\in{\mathcal{G}}\setminus\\{\textsf{group}(\ell)\\}$; then,
${\left\lceil\frac{|G\setminus
A(t)|}{k(G)}\right\rceil}={\left\lceil\frac{|G\setminus
A(s)|}{k(G)}\right\rceil}\leq p({\mathcal{I}})-1$. The equality is by Step 10
and the inequality is by the assumption on $A(s)$. Otherwise,

${\left\lceil\frac{|G\setminus
A(t)|}{k(G)}\right\rceil}={\left\lceil\frac{|G\setminus A(s)|+|\\{\ell\\}\cap
G|-|\\{y\\}\cap G|}{k(G)}\right\rceil}={\left\lceil\frac{|G\setminus
A(s)|}{k(G)}\right\rceil}\leq p({\mathcal{I}})-1$

The equalities are by Step 10 and the inequality is by the assumption on
$A(s)$. Therefore, by Claim I.2 and the above loop invariant the claim
follows.

∎

######  Claim I.4.

$V({\mathcal{I}}\setminus A)\leq p({\mathcal{I}})-1$.

###### Proof.

$V({\mathcal{I}}\setminus
A)=\max_{G\in{\mathcal{G}}_{A}}{\left\lceil\frac{|G|}{k_{A}(G)}\right\rceil}=\max_{G\in{\mathcal{G}}}{\left\lceil\frac{|G\setminus
A|}{k(G)}\right\rceil}\leq p({\mathcal{I}})-1$

The second equality is by (102). The last inequality is because for all
bounding groups $G\in{\mathcal{B}}({\mathcal{I}})$ it holds that
${\left\lceil\frac{|G\setminus A|}{k(G)}\right\rceil}\leq p({\mathcal{I}})-1$
by Claim I.3 and for all
$G\in{\mathcal{G}}\setminus{\mathcal{B}}({\mathcal{I}})$ it holds that
${\left\lceil\frac{|G\setminus A|}{k(G)}\right\rceil}\leq p({\mathcal{I}})-1$
by (104). ∎

######  Lemma I.5.

If $s(A)<1-\delta$, then Algorithm 10 returns a packing for ${\mathcal{I}}$
with at most $p({\mathcal{I}})$ bins.

###### Proof.

We use the following auxiliary claims.

######  Claim I.6.

For all $G\in{\mathcal{G}}$ it holds that $|A\cap G|=\min\\{k(G),|G|\\}$.

###### Proof.

Assume towards a contradiction that there is $G\in{\mathcal{G}}$ such that
$|A\cap G|<\min\\{k(G),|G|\\}$. Therefore, there is $\ell\in G\setminus A$.
Observe $A$ after Step 10. Because $s(A)<1-\delta$, it follows by Step 10 and
Step 10 that $\ell$ can be added to $A$ in contradiction. ∎

######  Claim I.7.

For all $G\in{\mathcal{G}}$, $\ell\in A\cap G$ and $y\in G\setminus A$ it
holds that $s(y)\leq s(\ell)$.

###### Proof.

Assume towards a contradiction that there are $\ell\in A,y\in I\setminus A$
such that $s(y)>s(\ell)$ and
$\textnormal{{group}}(\ell)=\textnormal{{group}}(y)$. Observe $A$ after Step
10. Because $s(A)<1-\delta$, it follows by Step 10 and Step 10 that $y$ can be
added to $A$ and $\ell$ can be removed from $A$ in contradiction.

∎

######  Claim I.8.

If $s(A)<1-\delta$ then $V({\mathcal{I}}\setminus A)\leq V({\mathcal{I}})-1$.

###### Proof.

$\displaystyle V({\mathcal{I}}\setminus A)={}$
$\displaystyle\max_{G\in{\mathcal{G}}_{A}}{\left\lceil\frac{|G|}{k_{A}(G)}\right\rceil}=\max_{G\in{\mathcal{G}}}{\left\lceil\frac{|G\setminus
A|}{k(G)}\right\rceil}=\max_{G\in{\mathcal{G}}}{\left\lceil\frac{|G|-\min\\{|G|,k(G)\\}}{k(G)}\right\rceil}$
$\displaystyle\leq{}$
$\displaystyle\max_{G\in{\mathcal{G}}}{\left\lceil\frac{|G|}{k(G)}-1\right\rceil}=\max_{G\in{\mathcal{G}}}{\left\lceil\frac{|G|}{k(G)}\right\rceil}-1=V({\mathcal{I}})-1.$

The second equality is by (102). The third equality is by Claim I.6. The
inequality is because if $\min\\{|G|,k(G)\\}=|G|$ then it follows that
$V({\mathcal{I}}\setminus A)=0\leq V({\mathcal{I}})-1$, assuming that
$I\neq\emptyset$. ∎

Now, back to the proof of Lemma I.5. We prove the lemma by induction on
$V({\mathcal{I}})$. For $V({\mathcal{I}})=1$, it follows that
$V({\mathcal{I}}\setminus A)\leq V({\mathcal{I}})-1$ by Claim I.8 and it
follows that $V({\mathcal{I}}\setminus A)=0$. Therefore, by (102) it follows
that $I\setminus A=\emptyset$ and therefore $A=I$. Thus, $I$ is a
configuration of ${\mathcal{I}}$ by Claim I.3. Hence, by Step 10 and Step 10
it holds that Algorithm 10 returns $(I)$, which is a packing of
${\mathcal{I}}$ because $I=A$ and $A$ is a configuration of ${\mathcal{I}}$.

Recall that for any BPP instance ${\mathcal{J}}$ we define $A^{{\mathcal{J}}}$
as the first entry of the tuple returned by
$\textsf{Greedy}({\mathcal{J}},\delta)$. Now, for the assumption of the
induction, assume that $V({\mathcal{I}})\geq 2$. Assume that for any BPP
instance ${\mathcal{J}}=(J,{\mathcal{G}}_{J},s_{J},k_{J})$ such that (i)
$s_{J}(\ell)\leq\delta~{}\forall\ell\in J$, (ii) $V({\mathcal{J}})\leq
V({\mathcal{I}})-1$, and (iii) $s\left(A^{{\mathcal{J}}}\right)<1-\delta$, it
holds that $\textsf{Greedy}({\mathcal{J}},\delta)$ returns a packing of
${\mathcal{J}}$ of at most $V({\mathcal{J}})$ bins. For the step of the
induction, we use the following auxiliary claim.

######  Claim I.9.

$s\left(A^{{\mathcal{I}}\setminus A}\right)<1-\delta$.

###### Proof.

Assume towards a contradiction that $s\left(A^{{\mathcal{I}}\setminus
A}\right)\geq 1-\delta$. Therefore, it holds that
$s\left(A^{{\mathcal{I}}\setminus A}\right)>s(A)$; thus, we reach a
contradiction if one of the following conditions hold.

  1. 1.

$|A^{{\mathcal{I}}\setminus A}|>|A|$. Therefore,

$|A^{{\mathcal{I}}\setminus
A}|=\sum_{G\in{\mathcal{G}}}|A^{{\mathcal{I}}\setminus A}\cap
G|\leq\sum_{G\in{\mathcal{G}}}\min\\{k(G),|G|\\}=\sum_{G\in{\mathcal{G}}}|A\cap
G|=|A|.$

The inequality is because $A^{{\mathcal{I}}\setminus
A}\in{\mathcal{C}}({\mathcal{I}}\setminus A)$ by Claim I.3 and
${\mathcal{C}}({\mathcal{I}}\setminus A)\subseteq{\mathcal{C}}({\mathcal{I}})$
by (102).999Claim I.3 is stated and proven for $A$ and not for
$A^{{\mathcal{I}}\setminus A}$ for simplicity; however, the exact arguments
can be used for proving the claim also for $A^{{\mathcal{I}}\setminus A}$. The
second equality is by Claim I.6. The above equation is a contradiction that
$|A^{{\mathcal{I}}\setminus A}|>|A|$.

  2. 2.

There are $G\in{\mathcal{G}}$, $\ell\in A^{{\mathcal{I}}\setminus A}\cap G$
and $y\in G\setminus A^{{\mathcal{I}}\setminus A}$ such that $s(y)>s(\ell)$.
This is a contradiction to Claim I.7.

Otherwise, if the two cases above are not satisfied:

$s\left(A^{{\mathcal{I}}\setminus
A}\right)=\sum_{G\in{\mathcal{G}}~{}}s\left(G\cap A^{{\mathcal{I}}\setminus
A}\right)\leq\sum_{G\in{\mathcal{G}}~{}}s\left(G\cap A\right)=s(A).$

The inequality is by Claim I.6, Claim I.7 and that conditions 1. and 2. above
are not satisfied. By the above equation, we reach a contradiction that
$s\left(A^{{\mathcal{I}}\setminus A}\right)>s(A)$.

∎

By Claim I.9 we can apply the assumption of the induction on
${\mathcal{I}}\setminus A$; thus, Algorithm 10 returns a packing for
${\mathcal{I}}\setminus A$ with at most $V({\mathcal{I}}\setminus A)$ bins. In
addition, by Claim I.8 it holds that $V({\mathcal{I}}\setminus A)\leq
V({\mathcal{I}})-1$. Hence, by Step 10, the returned object by the algorithm
is $X=(A)\oplus\textsf{Greedy}\left(\mathcal{I}\setminus A,\delta\right)$
which contains at most $V({\mathcal{I}})$ entries. Moreover, it holds that $X$
is a packing of ${\mathcal{I}}$ by Claim H.1 using the following arguments:
(i) $A\in{\mathcal{C}}({\mathcal{I}})$ by Claim I.3 and (ii)
$\textsf{Greedy}\left(\mathcal{I}\setminus A,\delta\right)$ is a packing of
${\mathcal{I}}\setminus A$ by the assumption of the induction. ∎

######  Lemma I.10.

If $s(A)\geq 1-\delta$, then Algorithm 10 returns a packing for
${\mathcal{I}}$ with at most $p({\mathcal{I}})$ bins.

###### Proof.

We prove the claim by induction on $p({\mathcal{I}})$. For
$p({\mathcal{I}})\leq 1$ it holds by (103) that $I$ is a configuration of
${\mathcal{I}}$: it follows that $V({\mathcal{I}}\setminus A)=0$ by Claim I.4
and therefore, by (102) it follows that $I\setminus A=\emptyset$ and $A=I$.
Hence, by Step 10 and Step 10 it holds that Algorithm 10 returns $(I)$, which
is a packing of ${\mathcal{I}}$ because $I$ is a configuration of
${\mathcal{I}}$. Now, for the assumption of the induction, assume that for any
BPP instance ${\mathcal{J}}=(J,{\mathcal{G}}_{J},s_{J},k_{J})$ such that (i)
$s(\ell)\leq\delta~{}\forall\ell\in J$, (ii)
$s\left(A^{{\mathcal{J}}}\right)\geq 1-\delta$, and (iii)
$p({\mathcal{J}})\leq p({\mathcal{I}})-1$, it holds that
$\textsf{Greedy}({\mathcal{J}},\delta)$ returns a packing of ${\mathcal{J}}$
of at most $p({\mathcal{J}})$ bins.

For the step of the induction, assume that $p({\mathcal{I}})>1$. We use the
following auxiliary claim.

######  Claim I.11.

$p\left({\mathcal{I}}\setminus A\right)\leq p({\mathcal{I}})-1$.

###### Proof.

We use the following inequality.

$(1+2\delta)\cdot s(I\setminus
A)+2\leq(1+2\delta)\cdot(s(I)-(1-\delta))+2\leq(1+2\delta)\cdot s(I)-1+2\leq
p({\mathcal{I}})-1.$ (108)

The first inequality is because $s(A)\geq 1-\delta$. The second inequality is
because $\delta\in(0,0.5)$. The last inequality is by (103). Now,

$p\left({\mathcal{I}}\setminus A\right)=\max\left\\{(1+2\delta)\cdot
s(I\setminus A)+2,V({\mathcal{I}}\setminus A)\right\\}\leq
p({\mathcal{I}})-1.$

The first equality is by (102) and (103). The first inequality is by (108) and
Claim I.4 ∎

To conclude the lemma, observe the following two complementary cases.

  1. 1.

$s\left(A^{{\mathcal{I}}\setminus A}\right)\geq 1-\delta$. Then, by Claim I.11
and that $s\left(A^{{\mathcal{I}}\setminus A}\right)\geq 1-\delta$, we can
apply the assumption of the induction on ${\mathcal{I}}\setminus A$; thus,
Algorithm 10 applied on ${\mathcal{I}}\setminus A$ returns a packing for
${\mathcal{I}}\setminus A$ with at most $p\left({\mathcal{I}}\setminus
A\right)$ bins. In addition, by Claim I.11 it holds that
$p\left({\mathcal{I}}\setminus A\right)\leq p({\mathcal{I}})-1$. Hence, by
Step 10, the returned object by the algorithm on ${\mathcal{I}}$ is
$X=(A)\oplus\textsf{Greedy}\left(\mathcal{I}\setminus A,\delta\right)$ which
contains at most $1+p({\mathcal{I}})-1=p({\mathcal{I}})$ entries. Moreover, it
holds that $X$ is a packing of ${\mathcal{I}}$ by Claim H.1 using the
following arguments: (i) $A\in{\mathcal{C}}({\mathcal{I}})$ by Claim I.3 and
(ii) $\textsf{Greedy}\left(\mathcal{I}\setminus A,\delta\right)$ is a packing
of ${\mathcal{I}}\setminus A$ by the assumption of the induction.

  2. 2.

$s\left(A^{{\mathcal{I}}\setminus A}\right)<1-\delta$. Then, by Lemma I.5 and
that $s\left(A^{{\mathcal{I}}\setminus A}\right)<1-\delta$, Algorithm 10
applied on ${\mathcal{I}}\setminus A$ returns a packing for
${\mathcal{I}}\setminus A$ with at most $V({\mathcal{I}}\setminus A)$ bins. In
addition, by Claim I.4 it holds that $V({\mathcal{I}}\setminus A)\leq
p({\mathcal{I}})-1$. Hence, by Step 10, the returned object by the algorithm
is $X=(A)\oplus\textsf{Greedy}\left(\mathcal{I}\setminus A,\delta\right)$
which contains at most $1+p({\mathcal{I}})-1\leq p({\mathcal{I}})$ entries by
(103). Moreover, it holds that $X$ is a packing of ${\mathcal{I}}$ by Claim
H.1 using the following arguments: (i) $A\in{\mathcal{C}}({\mathcal{I}})$ by
Claim I.3 and (ii) $\textsf{Greedy}\left(\mathcal{I}\setminus A,\delta\right)$
is a packing of ${\mathcal{I}}\setminus A$ by Lemma I.5.

∎

######  Lemma I.12.

The running time of Algorithm 10 is $\textnormal{poly}(|{\mathcal{I}}|)$.

###### Proof.

Step 10 and Step 10 can be trivially computed in linear time in
$|{\mathcal{I}}|$. Step 10 can be computed in polynomial time in
$|{\mathcal{I}}|$ by (i) finding the bounding groups by (104), and (ii)
compute a bounding subset of ${\mathcal{I}}$ as follows. For each bounding
group $G$ we compute a bounding subset of $G$ by sorting the items and
choosing a minimal number of items according to the sorted order that form a
bounding subset of $G$ (for more details see the proof of Lemma H.2).

the while loop of Step 10 runs at most $|I|$ times because in each iteration
an item is added to the constructed bin and is not added again by Step 10 and
Step 10. Finally, the while loop of Step 10 runs at most ${|I|\choose
2}=\textnormal{poly}(|{\mathcal{I}}|)$ times because in each iteration two
items are replaced in Step 10, one is add to the constructed bin and the other
one is removed; for any two items, this can happen at most once by Step 10. We
conclude that the running time of Algorithm 10 is
$\textnormal{poly}(|{\mathcal{I}}|)$. ∎

Proof of Lemma 1.2: If $s(A)<1-\delta$ it holds that
$X=(A)\oplus\textsf{Greedy}\left(\mathcal{I}\setminus A,\delta\right)$ is a
packing of ${\mathcal{I}}$ of at most $V({\mathcal{I}})$ bins by Lemma I.5;
therefore, the number of bins is at most $p({\mathcal{I}})$ by (103).
Otherwise, it holds that $s(A)\geq 1-\delta$; thus, it holds that $X$ is a
packing of ${\mathcal{I}}$ of at most $p({\mathcal{I}})$ bins by Lemma I.10.
Therefore, the number of bins in the packing is at most
$p({\mathcal{I}})\leq(1+2\delta)\cdot\max\left\\{s(I),~{}V(\mathcal{I})\right\\}+2$
by (103). Finally, the running time is $\textnormal{poly}(|{\mathcal{I}}|)$ by
Lemma I.12.∎

## References

  * [1] Adany, R., Feldman, M., Haramaty, E., Khandekar, R., Schieber, B., Schwartz, R., Shachnai, H., Tamir, T.: All-or-nothing generalized assignment with application to scheduling advertising campaigns. ACM Transactions on Algorithms (TALG) 12(3), 1–25 (2016)
  * [2] Alon, N., Azar, Y., Woeginger, G.J., Yadid, T.: Approximation schemes for scheduling on parallel machines. Journal of Scheduling 1(1), 55–66 (1998)
  * [3] Anderson, D.P.: Boinc: A system for public-resource computing and storage. In: Fifth IEEE/ACM international workshop on grid computing. pp. 4–10. IEEE (2004)
  * [4] Anderson, D.P.: BOINC: A platform for volunteer computing. Journal of Grid Computing 18, 99 – 122 (2017)
  * [5] Bansal, N., Caprara, A., Sviridenko, M.: A new approximation method for set covering problems, with applications to multidimensional bin packing. SIAM Journal on Computing 39(4), 1256–1278 (2010)
  * [6] Bansal, N., Eliáš, M., Khan, A.: Improved approximation for vector bin packing. In: Proceedings of the twenty-seventh annual ACM-SIAM symposium on discrete algorithms. pp. 1561–1579. SIAM (2016)
  * [7] Caprara, A., Kellerer, H., Pferschy, U.: Approximation schemes for ordered vector packing problems. Naval Research Logistics (NRL) 50(1), 58–69 (2003)
  * [8] Chekuri, C., Khanna, S.: A polynomial time approximation scheme for the multiple knapsack problem. SIAM Journal on Computing 35(3), 713–728 (2005)
  * [9] Christensen, H.I., Khan, A., Pokutta, S., Tetali, P.: Approximation and online algorithms for multidimensional bin packing: A survey. Computer Science Review 24, 63–79 (2017)
  * [10] Coffman, E.G., Csirik, J., Galambos, G., Martello, S., Vigo, D.: Bin packing approximation algorithms: survey and classification. In: Handbook of combinatorial optimization, pp. 455–531 (2013)
  * [11] Das, S., Wiese, A.: On minimizing the makespan when some jobs cannot be assigned on the same machine. In: 25th Annual European Symposium on Algorithms, ESA. pp. 31:1–31:14 (2017)
  * [12] Doron-Arad, I., Kulik, A., Shachnai, H.: An APTAS for bin packing with clique-graph conflicts. arXiv preprint arXiv:2011.04273 (2020)
  * [13] Doron-Arad, I., Kulik, A., Shachnai, H.: An APTAS for bin packing with clique-graph conflicts. In: Algorithms and Data Structures - 17th International Symposium, WADS 2021, Virtual Event, August 9-11, 2021, Proceedings. pp. 286–299 (2021)
  * [14] Doron-Arad, I., Kulik, A., Shachnai, H.: An fptas for budgeted laminar matroid independent set. arXiv preprint arXiv:2304.13984 (2023)
  * [15] Edmonds, J.: Minimum partition of a matroid into independent subsets. J. Res. Nat. Bur. Standards Sect. B 69, 67–72 (1965)
  * [16] Epstein, L., Levin, A.: AFPTAS results for common variants of bin packing: A new method for handling the small items. SIAM Journal on Optimization 20(6), 3121–3145 (2010)
  * [17] Fairstein, Y., Kulik, A., Shachnai, H.: Modular and submodular optimization with multiple knapsack constraints via fractional grouping. In: 29th Annual European Symposium on Algorithms, ESA 2021, September 6-8, 2021, Lisbon, Portugal (Virtual Conference). pp. 41:1–41:16 (2021)
  * [18] Fernandez de la Vega, W., Lueker, G.S.: Bin packing can be solved within 1 + $\varepsilon$ in linear time. Combinatorica 1, 349–355 (1981)
  * [19] Fleischer, L., Goemans, M.X., Mirrokni, V.S., Sviridenko, M.: Tight approximation algorithms for maximum separable assignment problems. Mathematics of Operations Research 36(3), 416–431 (2011)
  * [20] Gabow, H.N., Westermann, H.H.: Forests, frames, and games: algorithms for matroid sums and applications. Algorithmica 7(1), 465–497 (1992)
  * [21] Grage, K., Jansen, K., Klein, K.M.: An EPTAS for machine scheduling with bag-constraints. In: The 31st ACM Symposium on Parallelism in Algorithms and Architectures. pp. 135–144 (2019)
  * [22] Grigoriadis, M.D., Khachiyan, L.G., Porkolab, L., Villavicencio, J.: Approximate max-min resource sharing for structured concave optimization. SIAM Journal on Optimization 11(4), 1081–1091 (2001)
  * [23] Grötschel, M., Lovász, L., Schrijver, A.: Geometric algorithms and combinatorial optimization, vol. 2. Springer Science & Business Media (2012)
  * [24] Guerine, M., Stockinger, M.B., Rosseti, I., Simonetti, L.G., Ocaña, K.A., Plastino, A., de Oliveira, D.: A provenance-based heuristic for preserving results confidentiality in cloud-based scientific workflows. Future Generation Computer Systems 97, 697–713 (2019)
  * [25] Hoberg, R., Rothvoss, T.: A logarithmic additive integrality gap for bin packing. In: Proceedings of the Twenty-Eighth Annual ACM-SIAM Symposium on Discrete Algorithms. pp. 2616–2625. SIAM (2017)
  * [26] Hochbaum, D.S. (ed.): Approximation Algorithms for NP-Hard Problems. PWS Publishing Co., USA (1996)
  * [27] Hochbaum, D.S., Shmoys, D.B.: Using dual approximation algorithms for scheduling problems theoretical and practical results. Journal of the ACM 34(1), 144–162 (1987)
  * [28] Hoffman, A.J., Kruskal, J.B.: Integral boundary points of convex polyhedra. In: Linear Inequalities and Related Systems.(AM-38), Volume 38, pp. 223–246. Princeton University Press (1956)
  * [29] Jansen, K.: An approximation scheme for bin packing with conflicts. Journal of combinatorial optimization 3(4), 363–377 (1999)
  * [30] Jansen, K.: An EPTAS for scheduling jobs on uniform processors: using an MILP relaxation with a constant number of integral variables. SIAM Journal on Discrete Mathematics 24(2), 457–485 (2010)
  * [31] Jansen, K.: Parameterized approximation scheme for the multiple knapsack problem. SIAM Journal on Computing 39(4), 1392–1412 (2010)
  * [32] Jansen, K., Klein, K., Verschae, J.: Closing the gap for makespan scheduling via sparsification techniques. In: 43rd International Colloquium on Automata, Languages, and Programming (ICALP),. pp. 72:1–72:13 (2016)
  * [33] Jansen, K., Öhring, S.R.: Approximation algorithms for time constrained scheduling. Inf. Comput. 132(2), 85–108 (1997)
  * [34] Karmarkar, N., Karp, R.M.: An efficient approximation scheme for the one-dimensional bin-packing problem. In: 23rd Annual Symposium on Foundations of Computer Science. pp. 312–320. IEEE (1982)
  * [35] Kellerer, H., Pferschy, U.: Cardinality constrained bin-packing problems. Annals of Operations Research 92, 335–348 (1999)
  * [36] Krause, K.L., Shen, V.Y., Schwetman, H.D.: Analysis of several task-scheduling algorithms for a model of multiprogramming computer systems. Journal of the ACM (JACM) 22(4), 522–550 (1975)
  * [37] Krause, K., Shen, V.Y., Schwetman, H.D.: Errata:“analysis of several task-scheduling algorithms for a model of multiprogramming computer systems”. Journal of the ACM (JACM) 24(3), 527 (1977)
  * [38] Kulik, A., Mnich, M., Shachnai, H.: An asymptotic $\left(\frac{4}{3}+\varepsilon\right)$-approximation for the 2-dimensional vector bin packing problem (2022). https://doi.org/10.48550/ARXIV.2205.12828, https://arxiv.org/abs/2205.12828
  * [39] Fernandez de La Vega, W., Lueker, G.S.: Bin packing can be solved within 1+ $\varepsilon$ in linear time. Combinatorica 1(4), 349–355 (1981)
  * [40] Leung, J.Y.: Bin packing with restricted piece sizes. Information Processing Letters 31(3), 145–149 (1989)
  * [41] Lin, Y., Shen, H.: Eafr: An energy-efficient adaptive file replication system in data-intensive clusters. IEEE Transactions on Parallel and Distributed Systems 28(4), 1017–1030 (2017)
  * [42] McCloskey, B., Shankar, A.: Approaches to bin packing with clique-graph conflicts. Computer Science Division, University of California (2005)
  * [43] Plotkin, S.A., Shmoys, D.B., Tardos, É.: Fast approximation algorithms for fractional packing and covering problems. Mathematics of Operations Research 20(2), 257–301 (1995)
  * [44] Rothvoß, T.: Approximating bin packing within O(log OPT * log log OPT) bins. In: 54th Annual IEEE Symposium on Foundations of Computer Science. pp. 20–29. IEEE Computer Society (2013)
  * [45] Schrijver, A.: Combinatorial optimization: polyhedra and efficiency, vol. 24. Springer (2003)
  * [46] Simchi-Levi, D.: New worst-case results for the bin-packing problem. Naval Research Logistics (NRL) 41(4), 579–585 (1994)
  * [47] Williamson, D.P., Shmoys, D.B.: The Design of Approximation Algorithms. Cambridge University Press (2011)