# What’s math got to do with patterns in fish?

Blake Shirman Department of Mathematical Sciences, DePaul University, Chicago,
IL, USA Alexandria Volkening Department of Mathematics, Purdue University,
West Lafayette, IN, USA Correspondence<EMAIL_ADDRESS>
###### Abstract

When you think of fish, what comes to mind? Maybe you think of pet goldfish,
movie characters like Dory or Nemo, or trout in a local river. One of the
things that all of these fish have in common is patterns in their skin. Nemo
sports black and white stripes in his orange skin, and trout have spots. Even
goldfish have a pattern — it’s just plain gold (and kinda boring). Why do some
fish have stripes, others have spots, and others have plain patterns? It turns
out that this is a tricky question, so scientists need tools from several
subjects to answer it. In this paper, we use biology, math, and computer
coding to help figure out how fish get different skin patterns.

### Keywords

math $|$ model $|$ zebrafish $|$ biology $|$ computer science $|$ coding $|$
pattern $|$ animal skin

## What are fish patterns made of?

Many types of animals have skin patterns. Zebras have stripes and leopards
have spots, but fish display all sorts of patterns. When you look at Figure 1,
what patterns do you see? There are fish with stripes, spots, mazes, and plain
patterns [2]. Which pattern is your favorite? Zebrafish (Figure 1g) are good
for doing experiments with, so a lot of scientists like zebrafish best [3, 6].
As you might guess from the word “zebra” in their name, zebrafish usually have
stripes.

When we look at zebrafish from far away, we see black and yellow stripes. If
we put a zebrafish under a microscope, however, we can see that its patterns
are made up of tiny dots, or pigment cells. Pigment cells come in different
sizes and colors. Some cells are big and black (biologists call them
“melanophores”). Other cells are small and yellow (their fancy name is
“xanthophores”). There are more colors of pigment cells out there, and even
humans have pigment cells in their skin.

Figure 1:  Fish have all sorts of different skin patterns, and here are some
examples. Their patterns are made of pigment cells. (Images from [2, 5, 3, 6]
and Wikimedia Commons; see end of document for sources.)

## Cells follow rules to create stripes

Okay, so fish patterns contain pigment cells. But what’s really cool is that
these cells actually have the job of making the stripes [8]! In this
movie222Watch a fish pattern form here:
https://www.pnas.org/content/104/12/4790/tab-figures-data. Once you click the
website, scroll down to “Supporting Information” and then click “SI Movie 1.”
[8] of a fish growing, do you see the black cells moving around? There are
small yellow cells there too, but they are harder to see. You can think of
each pigment cell in zebrafish skin like it’s a person moving in a group. You
can move around in the room, and cells can move around in the skin, too.

Black and yellow cells move and act in specific ways (or follow rules) to
create patterns. Yellow cells following black cells is one of these rules.
It’s almost like they are playing _tag, you’re it_! You can see a black cell
leading the way and a yellow cell following it in this movie333Watch yellow
cells following black cells here: https://www.youtube.com/watch?v=0wECUnwgN8A
[9]. Crazy!

## How do scientists study fish patterns?

Figure 2:  Cells follow rules to create patterns. To simulate patterns with a
mathematical model on a computer, think of fish skin as a checkerboard [1].
Then pick a rule for what color of cell goes in each square. Our rule is “One
side different, the rest the same.” If you put your finger on a cell in the
center (indicated by the pink square), one of its neighbors (indicated by
green squares) should have a different color than it, and the rest should be
the same color as it.

While some of the rules for how cells behave in zebrafish are already known
(like yellow cells following black cells [9]), many rules are a mystery.
Scientists are working hard to discover the rules now. Cells are tiny, and
it’s tough to see them and catch them interacting with each other in the fish
skin. Zebrafish also take a long time to grow (a few months) [6], so
biologists have to be patient with their experiments.

This is where math and computer coding enter the picture. Math is all about
patterns, rules, and puzzles, and it can be combined with biology to help
figure out how cells behave [1, 7]. Mathematicians and biologists use
different tools to study the same problems. Bringing their ideas together
makes the mystery of how zebrafish get their stripes easier to solve.

While biologists study fish by looking into microscopes and doing experiments,
applied mathematicians write down rules for how cells might behave [1, 7].
These rules and equations are called a mathematical model— they are a way of
describing something in a mathematical way. It is important to remember that
models are a good guess, but they could be wrong. To test their rules,
mathematicians use computer coding to simulate fish patterns. This is like
growing fish on a computer! One good thing about this is that simulating fish
patterns only takes a few minutes, not a few months. (But we wouldn’t
recommend trying to eat a simulated fish for dinner! :-)

Table 1:  Determine the option that follows the mathematical model in Figure
2. Point your finger at the pink square and then use it and its green
neighbors to choose the color of the X-marked square. In problem 6, the pink
square only has 3 neighbors because it is at the checkerboard boundary.
Answers: 1b, 2b, 3a, 4c, 5c (one of the focus cell’s neighbors must be yellow,
but it could be its right or bottom side), 6a (this cell is at the edge of the
checkerboard and only has 3 sides, so we have enough information now), 7b, 8b,
9b, 10b.

## Building a mathematical model of zebrafish stripes

Now we are ready to build a mathematical model of zebrafish patterns! Our
first step is to describe fish skin as a checkerboard (like in Figure 2). Each
square in our checkerboard can either be empty (white), have a yellow cell, or
have a black cell.

To simulate fish patterns, we need to tell the black and yellow cells where to
appear in the checkerboard. In other words, we need to describe their rules of
behavior. Biologists know that black cells mostly like to be near other black
cells, and yellow cells like to have mostly yellow cells around them. We use
this fact to build our mathematical model.

In Figure 2, each cell in the checkerboard has 4 sides (unless the cell is at
a corner or an edge). Our rule is “One side different, the rest the same.”
This means that if you put your finger on a yellow cell, one of the cells that
it touches should be black (different) and the other cells should be yellow
(the same). You can see some examples of cells that follow this rule and cells
that don’t in Figure 2. To get some practice with this rule, try filling in
the worksheet in Table 1 (answers are at the bottom of the table).

Our model has only one rule, and the truth is that we are simplifying all the
ways that cells actually behave in zebrafish. But that’s okay! When
mathematicians make models, they start with simple rules and then build more
realistic ones. Simple models help scientists learn, and our model is a first
step. If you’d like to see more realistic models that mathematicians are
working on now, check out the scientific study [7].

## Our puzzle: Simulating zebrafish stripes

What pattern does the rule “One side different, the rest the same” create? You
can try this out and build a fish pattern using Figure 3. If you print Figure
3 and cut out the black and yellow cells, follow these steps:

  * •

Step 1: Zebrafish have one yellow stripe when they are young, so start by
putting 2 strips of yellow cells at the top of your fish checkerboard.

  * •

Step 2: Fill out the rest of the checkerboard by following the rule “One side
different, the rest the same.”

Start by filling out the squares at the top of the checkerboard and slowly
move your finger down the fish from left to right. What pattern forms?

Figure 3:  Use this image to create a fish pattern by combining biology and
math. Fill the checkerboard on the fish by following Steps 1 and 2. After you
finish, print this page again. It’s your turn to be the applied mathematician
now — make up your own rules to discover other patterns that you can create!
How would you change Step 2 to make a goldfish? You can also see how math is
combined with computer coding to simulate different patterns at this website.

By following the steps in Figure 3, you are doing the same steps by hand that
a computer would do to simulate a fish pattern. We’ve used computer coding to
create fish patterns — check out our website444Simulate a fish pattern here:
https://simulatingzebrafish.gitlab.io/zebrafish-outreach/ to simulate a fish
yourself!

## Why is studying fish patterns important, anyway?

Fish patterns might be nice to look at, but why do biologists and
mathematicians study them? That’s a really good question, and the answer has
to do with mutant patterns, medicine, and humans.

Turns out that zebrafish only _usually_ have stripes (hence the “zebra” part
of their name). As you know, stripes form because cells follow specific rules.
But what happens if cells don’t follow their normal rules? This can happen if
a fish has a mutation in its genes [3, 4]. In Figure 1, image (g) is not the
only picture of a zebrafish. Believe it or not, images (c) and (d) are also
zebrafish, but they are mutants!

Mutant zebrafish have different patterns because their cells follow different
rules. A fish’s genes specify these rules, so mutant genes lead to mutant
rules and mutant patterns. Humans can also have mutations in their genes, and
this can cause disease.

Zebrafish and humans look very different, but they actually have a lot of
similar genes. And that’s why studying zebrafish is important. If we can
figure out what rules cells follow to create normal and mutant patterns in
fish, we might be able to better understand normal and unhealthy behavior of
cells in humans down the line too.

## Now you get to be the applied mathematician!

Scientists are working hard to discover the ways that black and yellow cells
behave in normal zebrafish, and to figure out how these cells change their
behavior in mutant fish [4, 3]. For mathematicians who study zebrafish, their
job is to think of a fish pattern (like stripes) and then find the rules that
can create that pattern [1, 7]. This is a puzzle.

How would you change the rule “One side different, the rest the same” to
create a plain yellow fish? What rules and starting points create diagonal
stripes? Try out your own rules by printing Figure 3 out again. By working on
these questions, you are taking on the role of mathematician yourself! What
fish patterns can you create with math?

## Glossary

genes

    The biological instructions that an animal gets from their parents to specify their features (like stripe or spot patterns for a fish, or green or brown eyes for a person)
mathematical model

    Rules or equations that describe something in a mathematical way. In this paper, the rule “One side different, the rest the same” is part of our mathematical model
mutation

     A change in an animal’s genes from the normal setting to something else. Zebrafish usually have stripes (Figure 1g), but mutant zebrafish (Figures 1c, 1d) have maze or spot patterns
pigment cell

    Cells are small units that all living things are made of (we are made of cells, and so is fish skin). Pigment cells are special cells that contain pigment (color)
simulate

     To imitate a real-world system on a computer. You can simulate fish patterns on this website111Simulate a fish pattern here: https://simulatingzebrafish.gitlab.io/zebrafish-outreach/ with computer coding and a mathematical model, which provides directions to the computer

## References

  * [1] D. Bullara and Y. De Decker, Pigment cell movement is not required for generation of Turing patterns in zebrafish skin, Nat Commun, 6 (2015).
  * [2] H. G. Frohnhöfer, J. Krauss, H. M. Maischein, and C. Nüsslein-Volhard, Iridophores and their interactions with other chromatophores are required for stripe formation in zebrafish, Development, 140 (2013), pp. 2997–3007.
  * [3] U. Irion and C. Nüsslein-Volhard, The identification of genes involved in the evolution of color patterns in fish, Curr Opin Genet Dev, 57 (2019), pp. 31–38.
  * [4] L. B. Patterson and D. M. Parichy, Zebrafish pigment pattern formation: Insights into the development and evolution of adult form, Annu Rev Genet, 53 (2019), pp. 505–530.
  * [5] P. Salis, T. Lorin, V. Lewis, C. Rey, A. Marcionetti, M. L. Escande, N. Roux, L. Besseau, N. Salamin, M. Sémon, D. Parichy, J. N. Volff, and V. Laudet, Developmental and comparative transcriptomic identification of iridophore contribution to white barring in clownfish, Pigment Cell Melanoma Res, 32 (2019), pp. 391–402.
  * [6] A. P. Singh and C. Nüsslein-Volhard, Zebrafish stripes as a model for vertebrate colour pattern formation, Curr Biol, 25 (2015), pp. R81–R92.
  * [7] A. Volkening and B. Sandstede, Iridophores as a source of robustness in zebrafish stripes and variability in _Danio_ patterns, Nat Commun, 9 (2018).
  * [8] M. Yamaguchi, E. Yoshimoto, and S. Kondo, Pattern regulation in the stripe of zebrafish suggests an underlying dynamic and autonomous mechanism, Proc Natl Acad Sci USA, 104 (2007), pp. 4790–4793.
  * [9] H. Yamanaka and S. Kondo, In vitro analysis suggests that difference in cell movement during direct interaction can generate various pigment patterns in vivo, Proc Natl Acad Sci USA, 111 (2014), pp. 1867–1872.

### Image Sources

Figure 1: (a) from [5] with permission from John Wiley & Sons, Copyright
(2019) John Wiley & Sons A/S; (b) from [3] (licensed under CC BY-NC-ND 4.0);
(c) from [2] (CC-BY 3.0); (d) from [6] with permission from Elsevier,
Copyright (2015) Elsevier Ltd.; (e) from Zouavman Le Zouave (CC BY-SA) via
Wikimedia Commons; (f) from Bjwebb at English Wikipedia, Public domain, via
Wikimedia Commons; (g) adapted and combined from [2] (CC-BY 3.0) and Oregon
State University (CC BY-SA) via Wikimedia Commons. Figure 2: “Cells in fish
skin” image adapted from [2] (CC-BY 3.0).

### Conflict of Interest Statement

We declare no conflict of interest.

### Acknowledgments

We are grateful to Jithin George, Katelyn Joy Leisman, Niall Mangan, and Sasha
Shirman for helpful feedback on this manuscript, and BS thanks Sasha for
providing suggestions on JavaScript. We also thank Domenico Bullara and
Yannick De Decker for publishing their JavaScript code as part of [1], since
it helped inspire the code that we developed for our website. AV has been
supported in part by the National Science Foundation under grant no.
DMS-1764421 and by the Simons Foundation/SFARI under grant no. 597491-RWC.

Author Biographies

Blake Shirman is a Master’s student at DePaul University, and he hopes to one
day get a Ph.D. in pure math. He spends his time looking at logic puzzles and
games to see if math can make things easier, sometimes in fun and unexpected
ways! One of his favorite puzzles is trying to predict cell patterns (like how
to place checkers on a checkerboard with special rules about where they go).
It wasn’t until he met Alexandria that Blake learned diagrams like this were
being used to describe the vibrant stripes and spots on many of our undersea
friends!

Alexandria Volkening is an Assistant Professor in the Department of
Mathematics at Purdue University, and she got her Ph.D. in applied mathematics
from Brown University. She studies all sorts of stuff with math (zebrafish
patterns, of course, but also political elections, social media, and crowds of
people — turns out lots of things follow rules of behavior that can be
described mathematically!). Her favorite fish to simulate (so far) is
zebrafish, but her favorite fish for dinner is salmon. And no, she doesn’t
have any pet zebrafish, so she always enjoys visiting biologists’ labs to see
them!


###### Abstract

We give a column generation based branch and bound algorithm for coalition
structure generation over graphs problem using valuation functions for which
this problem is proven to be NP-complete. For a given graph $G=(V,E)$ and a
valuation function $w:2^{V}\to{\mathbb{R}}$, the problem is to find the most
valuable coalition structure (or partition) of $V$. We consider two cases:
first when the value of a coalition is the sum of the weights of its edges
which can be positive or negative, second when the value of a coalition takes
account of both inter- and intra-coalitional disagreements and agreements,
respectively. For both valuations we give experimental results which cover for
the first time sets of more than forty agents.

For another valuation function (coordination) we give only the theoretical
considerations in the appendix.

Keywords: $coalition\;structure\;generation$, $linear\;programming$,
$quadratic\;programming$, $column\;generation$,
$coalition\;valuation\;function$.

††The publisher does not claim any copyright for the technical reports. The
author keeps the full copyright for the paper, and is thus free to transfer
the copyright to a publisher if the paper is accepted for publication
elsewhere.

Faculty of Computer Science, Alexandru Ioan Cuza University Iaşi

Technical Report TR 20-03, April 2020

A Branch and Bound Algorithm for Coalition Structure Generation over Graphs

Emanuel Florentin OLARIU, Cristian FRĂSINARU, Albert Abel POLICIUC
Faculty of Computer Science, Alexandru Ioan Cuza University Iaşi,

General Berthelot 16, 700483 Iasi̧, Romania,

Email<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
##  1 Introduction

Coalition structure generation (CSG) is a major problem in artificial
intelligence ([13]), multi-agent systems ([2], [14]) communication networks,
cooperative game theory ([3], [4], [12]), scheduling ([7]), economic theory
(like combinatorial auctions) etc. Given a set of agents
$V=\\{1,2,\ldots,n\\}$, and a valuation function $v:2^{V}\to{\mathbb{R}}$
assigning a value to any coalition of agents, the problem is to partition the
set of agents into pairwise disjoint classes (coalitions) such that the sum of
their values is maximized. CSG is one of the steps in the coalition formation
process [11] which may include also optimizing the coalitions performance and
rewarding the coalitions value among the members.

CSG comes from real-world applications: consider a set of agents who can
cooperate by working in coalitions. Some of them work better together while
others find difficult to cooperate. The problem is to maximize the so called
social well-fare or the total value of the designated coalitions. Classical
cooperative game theory usually uses a super-additive valuation function that
values better a merged coalition than the sum of of the values of the
component coalitions. This leads to the great coalition formation which can be
costly to coordinate and/or manipulate. Besides these considerations there are
natural constraints on possible coalitions, hence the abstractions like super-
additivity are not always appropriate for modeling the coalitions values and
the agents should be divided into smaller coalitions.

The problem is computationally difficult: the input specification for all
$2^{|V|}$ possible coalitions is intractable even for reasonable values of
$|V|$ and the computational difficulty maintains even under quite restrictive
assumptions. The literature presents a very large number of approaches to CSG
problem like (see [10]): dynamic programming, meta heuristic methods, branch
and bound algorithms based on dividing the searching space, anytime algorithms
that maintain a monotonically improving feasible solution etc.

Our setting is that of (mixed) integer linear/quadratic programming and is
based on the classical model of the set partitioning problem in a graph. Given
a graph $G=(V,E)$ the vertices are agents and the edges represent connections
between agents. We study two types of valuation functions: the edge sum and
the correlation functions. The value of a coalition by the edge sum valuation
function ([3]) is the sum of the weights of its edges (suppose we have a
weight function defined on the set of edges); this function frequently occurs
in communication networks and cooperative game theory. The correlation
valuation function ([1]) is defined on the entire coalition structure,
although it can be defined on coalitions only, and occurs in the clustering
framework. This function takes account of the agreements from inside and the
disagreements from outside the structure - an edge being labeled with a plus
or a minus depending on whether the involved agents are similar or different.

Our approach is based on solving an integer linear programming problem, which
is equivalent with CSG, using a branch and bound algorithm in which the
problems in nodes are built by means of column generation method. This
approach works as long as the involved sub-problem can be solved. We proved
that this can be done for edge sum and correlation valuation functions since
the sub-problem becomes a quadratic knapsack problem with forbidden
configurations. This variant of quadratic knapsack problem can be solved in
practice using the corresponding integer quadratic programming model or the
mixed integer linear programming equivalent model.

Our numerical results show that we can approach by this method sets of up to
$45$ agents depending on the magnitude of the weight function for edge sum
valuation function and up to $40$ agents for the correlation valuation
function, using a regular home PC.

The paper is organized as follows: section 2 describes the LP model and the
column generation framework, section 3 is dedicated to the two valuation
functions outlining the corresponding sub-problems and their linearizations,
section 4 describes the branch and bound algorithm, and the last section
contains the numerical results and the conclusions. The appendix contains ..

##  2 LP model and column generation

Consider $V=\\{1,2,\ldots,n\\}$ to be a set of $n$ agents and
$v:2^{V}\to{\mathbb{R}}$ a valuation function on the power set of $V$. A
coalition structure is a collection of disjoint exhaustive subsets
$C_{1},C_{2},\ldots,C_{p}$, where $C_{i}\cap C_{j}=\varnothing$ for all $1\leq
i<j\leq p$ and $\displaystyle\bigcup_{i=1}^{p}C_{i}=V$.

The problem of finding a coalition structure of maximum value is equivalent
with the set partitioning problem (SPP) (see [7]):

$\displaystyle\begin{array}[]{rcl}&max&\displaystyle\sum_{j=0}^{2^{n}-1}v(C_{j})x_{j}\end{array}$
(2)
$\displaystyle\begin{array}[]{rcl}\displaystyle\sum_{j=0}^{2^{n}-1}a_{ij}x_{j}&=&1,\forall
i\in V,\end{array}$ (4) $\displaystyle\begin{array}[]{rcl}\displaystyle
x_{j}&\in&\\{0,1\\},\forall j\in\\{0,1,\ldots,2^{n}-1\\}\end{array}$ (6)

where $\\{C_{0},C_{1},\ldots,C_{2^{n}-1}\\}$ is an enumeration of $2^{V}$ and
$(a_{ij})_{i\in V}$ is the characteristic vector of $C_{j}$, for each $j$,
that is

$\displaystyle a_{ij}=\left\\{\begin{array}[]{rl}1,&\mbox{if }i\in C_{j}\\\
0,&\mbox{otherwise}\end{array}\right.,\mbox{ for each
}j\in\\{0,1,\ldots,2^{n}-1\\}.$ (9)

We can relax the integrality constraints by replacing (6) with

$\displaystyle\begin{array}[]{rcl}\displaystyle x_{j}&\geq&0,\forall
j\in\\{0,1,\ldots,2^{n}-1\\}\end{array}$ (11)

Writing the original relaxed problem as a minimum one means to replace (2) by
(13) (by ignoring the minus in front of min)

$\displaystyle\begin{array}[]{rcl}&min&\displaystyle\left(\sum_{j=0}^{2^{n}-1}-v(C_{j})x_{j}\right)\end{array}$
(13)

The dual of the problem (13), (4) , (11) is

$\displaystyle\begin{array}[]{rcl}&max&\displaystyle\left(\sum_{j=1}^{n}\pi_{i}\right)\end{array}$
(15)
$\displaystyle\begin{array}[]{rcl}\displaystyle\sum_{i=1}^{n}a_{ij}\pi_{i}&\leq&-v(C_{j}),\forall
j\in\\{0,1,\ldots,2^{n}-1\\}\end{array}$ (17)
$\displaystyle\begin{array}[]{rcl}\displaystyle\pi_{i}\in{\mathbb{R}},\forall
i\in V\end{array}$ (19)

Let $C_{0}=\varnothing$, and $C_{j}=\\{j\\}$, for any
$j\in\\{1,2,\ldots,m\\}$. An initial feasible basic solution to the problem
(13), (4), (11) could be that corresponding to the coalition structure
$\displaystyle\left(\\{1\\},\\{2\\},\ldots\\{m\\}\right)$, that is,
$x_{1},x_{2},\ldots,x_{m}$. Now we restrict the problem to a small number of
variables including $x_{1},x_{2},\ldots,x_{m}$ and we get the restricted
master problem (RMP):

$\displaystyle\begin{array}[]{rcl}&min&\displaystyle\left(\sum_{j\in
J}-v(C_{j})x_{j}\right)\end{array}$ (21)
$\displaystyle\begin{array}[]{rcl}\displaystyle\sum_{j\in
J}a_{ij}x_{j}&=&1,\forall i\in V,\end{array}$ (23)
$\displaystyle\begin{array}[]{rcl}\displaystyle x_{j}&\geq&0,\forall j\in
J\end{array}$ (25)

where $\displaystyle J=\left\\{1,2,\ldots,m\right\\}\cup
J^{\prime}\subseteq\\{0,1,\ldots,2^{m}-1\\}$.

The dual of the RMP is

$\displaystyle\begin{array}[]{rcl}&max&\displaystyle\left(\sum_{i=1}^{n}\pi_{i}\right)\end{array}$
(27)
$\displaystyle\begin{array}[]{rcl}\displaystyle\sum_{i=1}^{n}a_{ij}\pi_{i}&\leq&-v(C_{j}),\forall
j\in J\end{array}$ (29)
$\displaystyle\begin{array}[]{rcl}\displaystyle\pi_{i}\in{\mathbb{R}},\forall
1\leq i\leq n\end{array}$ (31)

Let $({\mathbf{x}},\mathbold{\pi})$ be an optimum primal-dual solution for
this pair of problems; we look for a non-basic variable (column) with the
minimum negative reduced cost (Dantzig rule) - by solving a corresponding sub-
problem - that would be added to the current restricted master problem. When
such a variable doesn’t exist we can stop: we have an optimum solution to the
primal problem. The sub-problem is to find

$\displaystyle\begin{array}[]{rcl}\displaystyle
j_{0}=\operatorname*{arg\,max}_{j\in\\{0,1,\ldots,2^{n}-1\\}\setminus
J}\left(v(C_{j})+\sum_{i=1}^{n}a_{ij}\pi_{i}\right)\end{array}$ (33)

The arising question is: how can we solve problem (33)? The answer depends on
the form of the valuation function $v(\cdot)$. In the following sections we
will analyze this sub-problem for two of the most frequent used valuation
functions defined for coalition structure over graphs: the edge-sum and the
correlation valuation functions.

##  3 Coalition Structure Generation over Graphs

###  3.1 The edge-sum valuation

Let $G=(V,E)$ be graph and $w:E\to{\mathbb{R}}$ a weight function on its
edges. The corresponding edge-sum coalition valuation function is

$v:{\mathcal{P}}(V)\to{\mathbb{R}},v(C)=\sum_{ij\in E:i,j\in C}w_{ij},\forall
C\subseteq V.$

This function was extensively studied in the context of cooperative game
theory ([3]) and the corresponding CSG problem is NP-hard being as hard as the
MAX-CUT problem.

Now, if ${\mathbf{v}}$ is the characteristic vector of the generic coalition
$C$, then the sub-problem becomes:

$\displaystyle\max_{{\mathbf{v}}\in\\{0,1\\}^{m},v\not=\chi_{C_{j}},\forall
j\in J}{\left(\sum_{ij\in E}w_{ij}v_{i}v_{j}+\sum_{i\in
V}\pi_{i}v_{i}\right)}>0$ (34)

That is, the sub-problem is a quadratic knapsack problem with forbidden
configurations (QKPf). We can find a new column to add to the restricted
master problem if and only if this quadratic knapsack problem has a strictly
positive optimal objective value.

QKPf could be a computationally difficult problem since the quadratic knapsack
problem (QKP) is known to be NP-hard (being a generalization of Clique
problem). By rephrasing it we get:

$\displaystyle\begin{array}[]{rcl}&max&\displaystyle\left(\sum_{ij\in
E}w_{ij}v_{i}v_{j}+\sum_{i\in V}\pi_{i}v_{i}\right)\end{array}$ (36)
$\displaystyle\begin{array}[]{rcl}\displaystyle\sum_{i\in
C_{j}}v_{i}+\sum_{i\notin C_{j}}(1-v_{i})&\geq&1,\forall j\in J\end{array}$
(38) $\displaystyle\begin{array}[]{rcl}\displaystyle
v_{i}&\in&\\{0,1\\},\forall i\in V\end{array}$ (40)

This problem has a quadratic objective but only linear constraints: an Integer
Quadratic Programming (IQP) problem. Such a problem can be solved by using a
mathematical optimization solver as such or by linearizing it first in order
to transform it into a Mixed Integer Linear Programming (MILP) problem. We
note here that the known methods used for solving such a problem are not
applicable here since the usual requirements are to have natural or positive
coefficients for the quadratic (some times non-diagonal) terms ([5], [8]).
Hence we had to settle to solve it as a general IQP problem using a
mathematical optimization solver.

###  3.2 The correlation valuation

Suppose we have a function $e:E\to\\{+,-\\}$ that assigns to each edge $ij\in
E$ the labels $+$ or $-$ depending on whether agents $i$ and $j$ are similar
or different (this function arises in clustering frameworks). A valuation
function that takes account of both inter- and intra-coalitional similarities
can be defined in the following way ([1]): we want to minimize the number of
mistakes: a positive mistake occurs when $e(ij)=-$, with $i$ and $j$ belonging
to the same coalition, a negative one occurs when $e(ij)=+$, but $i$ and $j$
belong to different coalitions. Or, equivalently, we maximize the number of
agreements; first define, for a coalition $C\subseteq V$, the intra- and
inter-coalitional connections

$Intra^{+}(C)=|\\{ij\in E\>:\>e(ij)=+,i,j\in C\\}|,$ $Inter^{-}(C)=|\\{ij\in
E\>:\>e(ij)=-,i\in C,j\notin C\\}\\}|.$

The correlation valuation function is

$v:{\mathcal{P}}(V)\to{\mathbb{R}},v(C)=Intra^{+}(C)+Inter^{-}(C)/2,\forall
C\subseteq V.$

It was proven in ([1]) that the corresponding sub-problem is $NP$-complete. We
define two weight functions on the set of edges of $G$,
$w^{+},w^{-}:E\to\\{0,1\\}$, by

$w^{+}(ij)=w^{+}_{ij}=\left\\{\begin{array}[]{rl}1,&\mbox{if }e(ij)=+\\\
0,&\mbox{if }e(ij)=-\end{array}\right.,$
$w^{-}(ij)=w^{-}_{ij}=\left\\{\begin{array}[]{rl}1,&\mbox{if }e(ij)=-\\\
0,&\mbox{if }e(ij)=+\end{array}\right.,\forall ij\in E.$

Let ${\mathbf{v}}$ be the characteristic vector of the generic coalition $C$,
then

$\displaystyle\begin{array}[]{l}\displaystyle v(C)=\frac{1}{2}\sum_{i\in
V}\sum_{j\in V}w^{+}_{ij}v_{i}v_{j}+\frac{1}{2}\sum_{i\in V}\sum_{j\in
V}w^{-}_{ij}v_{i}(1-v_{j}).\end{array}$ (42)

The sub-problem (33) becomes the following QKPf problem:

$\displaystyle\begin{array}[]{rcl}&max&\displaystyle\left(\frac{1}{2}\sum_{i\in
V}\sum_{j\in V}w^{+}_{ij}v_{i}v_{j}+\frac{1}{2}\sum_{i\in V}\sum_{j\in
V}w^{-}_{ij}v_{i}(1-v_{j})+\sum_{i\in V}\pi_{i}v_{i}\right)\end{array}$ (44)
$\displaystyle\begin{array}[]{rcl}\displaystyle\sum_{i\in
C_{j}}v_{i}+\sum_{i\notin C_{j}}(1-v_{i})&\geq&1,\forall j\in J\end{array}$
(46) $\displaystyle\begin{array}[]{rcl}\displaystyle
v_{i}&\in&\\{0,1\\},\forall i\in V\end{array}$ (48)

###  3.3 The linearization of QKPf

We present here the linearization due to Glover ([6]) for solving QKP. First,
we consider the following equivalent form of the objective function quadratic
fragment $\displaystyle\sum_{i=1}^{n}\sum_{j=1}^{n}w^{\prime}_{ij}v_{i}v_{j}$,
where $2w^{\prime}_{ij}=w_{ij}$, for each $i\not=j$, and $w^{\prime}_{ii}=0$,
for all $i\in V$.
Second,$\displaystyle\sum_{j=1}^{n}w^{\prime}_{ij}v_{i}v_{j}$ is replaced by a
new real variable $u_{i}$, and we add four new constraints

$\displaystyle\begin{array}[]{l}W_{i}^{{}^{\prime}-}v_{i}\leq u_{i}\leq
W_{i}^{{}^{\prime}+}v_{i},\\\
\displaystyle\sum_{j=1}^{n}w^{\prime}_{ij}v_{j}-W_{i}^{{}^{\prime}+}(1-v_{i})\leq
u_{i}\leq\displaystyle\sum_{j=1}^{n}w^{\prime}_{ij}v_{j}-W_{i}^{{}^{\prime}-}(1-v_{i}),\end{array}$
(51)

where $W_{i}^{{}^{\prime}-}$ and $W_{i}^{{}^{\prime}+}$ are a lower and,
respectively, an upper bound for
$\displaystyle\sum_{i=1}^{n}\sum_{j=1}^{n}w^{\prime}_{ij}v_{j}$, e.g.:

$\displaystyle\begin{array}[]{l}W_{i}^{{}^{\prime}-}=\displaystyle\sum_{j:w^{\prime}_{ij}<0}w^{\prime}_{ij}\mbox{
and
}W_{i}^{{}^{\prime}+}=\displaystyle\sum_{j:w^{\prime}_{ij}>0}w^{\prime}_{ij}\end{array}$
(53)

The problem (36) - (40) becomes, after replacing $w^{\prime}_{ij}$,
$W_{i}^{{}^{\prime}-}$, and $W_{i}^{{}^{\prime}+}$, by correspondingly
$w_{ij}$, $W_{i}^{-}$, and $W_{i}^{+}$ (and substituting $2u_{i}$ by $u_{i}$)

$\displaystyle\begin{array}[]{rcl}&max&\displaystyle\left(\sum_{i\in
V}\pi_{i}v_{i}+\sum_{i\in V}u_{i}\right)\end{array}$ (55)
$\displaystyle\begin{array}[]{rcl}\displaystyle\sum_{i\in
C_{j}}v_{i}+\sum_{i\notin C_{j}}(1-v_{i})&\geq&1,\forall j\in J,\end{array}$
(57) $\displaystyle\begin{array}[]{rcl}\displaystyle\sum_{j:ij\in
E}w_{ij}v_{j}+W_{i}^{+}v_{i}-u_{i}&\leq&W_{i}^{+},\forall i\in V,\end{array}$
(59) $\displaystyle\begin{array}[]{rcl}\displaystyle\sum_{j:ij\in
E}w_{ij}v_{j}+W_{i}^{-}v_{i}-u_{i}&\geq&W_{i}^{-},\forall i\in V,\end{array}$
(61) $\displaystyle\begin{array}[]{rcl}W_{i}^{-}v_{i}-u_{i}\leq 0,\forall i\in
V,\end{array}$ (63) $\displaystyle\begin{array}[]{rcl}W_{i}^{+}v_{i}-u_{i}\geq
0,\forall i\in V,\end{array}$ (65)
$\displaystyle\begin{array}[]{rcl}\displaystyle v_{i}&\in&\\{0,1\\},\forall
i\in V,\end{array}$ (67) $\displaystyle\begin{array}[]{rcl}\displaystyle
u_{i}&\in&{\mathbb{R}},\forall i\in V.\end{array}$ (69)

For the correlation valuation function the corresponding model can be
linearized as in a similar way; replacing the coefficients of the linear
fragment with $\displaystyle p_{i}=\frac{1}{2}\sum_{j\in
V}w^{-}_{ij}+\pi_{i},\forall i\in V$, we get the following linear formulation
of the problem (44) - (48):

$\displaystyle\begin{array}[]{rcl}&max&\displaystyle\left(\sum_{i\in
V}p_{i}v_{i}+\sum_{i\in V}u_{i}\right)\end{array}$ (71)
$\displaystyle\begin{array}[]{rcl}\displaystyle\sum_{i\in
C_{j}}v_{i}+\sum_{i\notin C_{j}}(1-v_{i})&\geq&1,\forall j\in J,\end{array}$
(73) $\displaystyle\begin{array}[]{rcl}\displaystyle\sum_{j:ij\in
E}w^{\prime}_{ij}v_{j}+W_{i}^{{}^{\prime}+}v_{i}-u_{i}&\leq&W_{i}^{{}^{\prime}+},\forall
i\in V,\end{array}$ (75)
$\displaystyle\begin{array}[]{rcl}\displaystyle\sum_{j:ij\in
E}w^{\prime}_{ij}v_{j}+W_{i}^{{}^{\prime}-}v_{i}-u_{i}&\geq&W_{i}^{{}^{\prime}-},\forall
i\in V,\end{array}$ (77)
$\displaystyle\begin{array}[]{rcl}W_{i}^{{}^{\prime}-}v_{i}-u_{i}\leq
0,\forall i\in V,\end{array}$ (79)
$\displaystyle\begin{array}[]{rcl}W_{i}^{{}^{\prime}+}v_{i}-u_{i}\geq
0,\forall i\in V,\end{array}$ (81)
$\displaystyle\begin{array}[]{rcl}\displaystyle v_{i}&\in&\\{0,1\\},\forall
i\in V,\end{array}$ (83) $\displaystyle\begin{array}[]{rcl}\displaystyle
u_{i}&\in&{\mathbb{R}},\forall i\in V,\end{array}$ (85)

where $w^{\prime}_{ij}=0.5\cdot(w^{+}_{ij}-w^{-}_{ij})$, for all $i,j\in V$,
and the upper and lower bounds ($W_{i}^{{}^{\prime}+}$ and
$W_{i}^{{}^{\prime}-}$) may have the same expressions like in (53).

These last problems are Mixed Integer Linear Programming (MILP) problems that
have a reasonable size and can theoretically be solved by the means of a
mathematical optimization solver. We presented these linear models for the
sake of completeness, but the numerical results for them are, by now, - in
terms of running time performance - not encouraging. The optimization solver
we used performs better on the quadratic models.

##  4 Branch and Bound algorithm

Suppose that the current node in the branching tree has a subset of already
covered agents, $\displaystyle U=\displaystyle\bigcup_{j\in J:x_{j}=1}C_{j}$,
and $U^{\prime}=V\setminus U$. For the edge-sum valuation the sub-problem
becomes:

$\displaystyle\begin{array}[]{rcl}&max&\displaystyle\left(\sum_{i\in
U^{\prime}}\sum_{j\in U^{\prime}}w_{ij}v_{i}v_{j}+\sum_{i\in
U^{\prime}}\pi_{i}v_{i}\right)\end{array}$ (87)
$\displaystyle\begin{array}[]{rcl}\displaystyle\sum_{i\in
C_{j}}v_{i}+\sum_{i\notin C_{j}}(1-v_{i})&\geq&1,\forall j\in J,\mbox{ s. t.
}C_{j}\subseteq U^{\prime}.\end{array}$ (89)
$\displaystyle\begin{array}[]{rcl}\displaystyle v_{i}&\in&\\{0,1\\},\forall
i\in U^{\prime}.\end{array}$ (91)

For the correlation valuation function the sub-problem (33) becomes:

$\displaystyle\begin{array}[]{rcl}&max&\displaystyle\left[\frac{1}{2}\sum_{i\in
U^{\prime}}\sum_{j\in U^{\prime}}(w^{+}_{ij}-w_{ij}^{-})v_{i}v_{j}+\sum_{i\in
U^{\prime}}\left(\frac{1}{2}\sum_{j\in
V}w^{-}_{ij}+\pi_{i}\right)v_{i}\right]\end{array}$ (93)
$\displaystyle\begin{array}[]{rcl}\displaystyle\sum_{i\in
C_{j}}v_{i}+\sum_{i\notin C_{j}}(1-v_{i})&\geq&1,\forall j\in J,\mbox{ s. t.
}C_{j}\subseteq U^{\prime}.\end{array}$ (95)
$\displaystyle\begin{array}[]{rcl}\displaystyle v_{i}&\in&\\{0,1\\},\forall
i\in U^{\prime}.\end{array}$ (97)

Both these problems have linearization variants.

At each node of the branching tree we first build the current LP relaxation in
two steps: (1) we reduce the variables number by taking account of the
branching variables along the path to the root, and (2) we add the necessary
variables to the parent node LP relaxation by the means of column generation
method.

Step (1) is implemented by effectively fixing to $0$ the variables $x_{h}$
such that $C_{h}\cap U\not=\varnothing$ (or, equivalently, $C_{h}\cap
C_{j}\not=\varnothing$, for some branching variable $x_{j}$ set to $1$), and
removing the branching variables $x_{i}$ set to 0. The fixing procedure can be
achieved by looking in each equation (4) that contains a branching variable
set to $1$. In this way any branching variable set to $1$, that corresponds to
a medium sized coalition, has the effect of drastically reducing the size of
the corresponding mathematical programming model.

While step (1) is basically the same for both valuation functions, step (2)
depends on the specific sub-problem. Step (2) consists in repeatedly solving
the corresponding sub-problem while the optimum objective function value is
strictly positive. The implementation, however, requires this value to be
positive within some tolerance. If this step would not have been subjected to
numerical restriction, then the algorithm would have been an exact method.

Our branching rule works in a classical way: we choose a variable, $x_{j}$,
from the optimal solution in the current node of the branching tree such that
its value is around $0.5$ (e. g. $x_{j}\in(0.35,0.65)$ \- if possible). It is
worth-noting that, for our specific problems, this enumerative method (branch
and bound) doesn’t need upper bounds because in almost all cases a feasible
solution occurs very early - mostly after performing one of the first type (1)
steps. The overall effect is that the branching tree has a medium size: we
limited the number of branching tree size to $40$, but this bound was hardly
reached.

##  5 Numerical results

In this section we evaluate the performance of our algorithm. A major part of
the running time of our algorithm is concentrated in the root node of the
branching tree, where step (2) adds hundreds of new variables to the root
model, hence hundreds of sub-problems to solve (fortunately Gurobi solver
quickly finds optimal solutions to these QKPf problems), while for the other
nodes of the tree the number of such sub-problems drastically reduces.

The algorithm has been written in Java and run on an Intel(R) Core (TM)
i5-7500 CPU 3.40GHz computer with 8GB RAM, under Ubuntu 18.04.4 LTS. The
linear and quadratic programming problems were solved using Gurobi 9.0 under
an Academic License.

Since there are no benchmarks in the literature for edge-sum or correlation
valuation functions, the algorithm has been tested on randomly generated
problems. Our test problems were built using the Gilbert model, that is, each
edge has a fixed probability of being present in the graph, independently of
any other edges (probability $1$ gives a complete graph). The weights on edges
are independently generated using a Gaussian distribution $N(0,0.2)$ (a
certain similar valuation function occurs in [9]). The name of the benchmark
file indicates the probability that an edge belongs to the graph (”p”), the
number of agents (”n”), and the number of the instance (”s”).

Table 1: Solutions and computational times for different instances. |  Instance | Solution | CPU time (seconds) | Number of | Gap
---|---|---|---|---|---
|  |  LP |  ILP |  overall |  root |  per node |  nodes |  variables |  int. sol.
Edge-sum valuation |  p0.8n40s0 |  17.877 |  17.866 |  775 |  547 |  258.3 |  3 |  316 |  1 | 0.06%
|  p0.8n40s1 |  21.495 |  21.353 |  1,833 |  1,027 |  122.2 |  15 |  849 |  1 | 0.66%
|  p0.8n40s2 |  20.108 |  20.108 |  720 |  720 |  720.0 |  1 |  253 |  1 | 0.00%
|  p0.8n40s3 |  19.327 |  19.290 |  812 |  650 |  162.4 |  5 |  320 |  1 | 0.19%
|  p0.8n40s4 |  18.166 |  18.099 |  897 |  584 |  179.4 |  5 |  318 |  1 | 0.36%
|  p0.8n45s0 |  22.606 |  22.411 |  6,672 |  3,565 |  166.8 |  40 |  1,598 |  5 | 0.86%
|  p0.8n45s1 |  22.580 |  22.480 |  4,451 |  2,902 |  404.63 |  11 |  622 |  2 | 0.44%
|  p0.8n45s2 |  21.485 |  21.432 |  2,943 |  1,805 |  196.2 |  15 |  799 |  2 | 0.24%
|  p0.8n45s3 |  23.884 |  23.884 |  2,815 |  2,815 |  2,815.0 |  1 |  312 |  1 | 0.00%
|  p0.8n45s4 |  21.500 |  21.500 |  4,472 |  4,472 |  4,472.0 |  1 |  284 |  1 | 0.00%
|  p1.0n35s0 |  17.154 |  17.154 |  346 |  346 |  346.0 |  1 |  203 |  1 | 0.00%
|  p1.0n35s1 |  16.427 |  16.383 |  446 |  344 |  89.2 |  5 |  273 |  1 | 0.26%
|  p1.0n35s2 |  19.406 |  19.241 |  1,651 |  590 |  47.1 |  35 |  1,231 |  2 | 0.85%
|  p1.0n35s3 |  15.500 |  15.388 |  630 |  398 |  90.0 |  7 |  334 |  2 | 0.72%
|  p1.0n35s4 |  19.064 |  18.987 |  979 |  414 |  108.7 |  9 |  515 |  1 | 0.40%
|  p1.0n40s0 |  25.286 |  25.286 |  3,860 |  3,860 |  3,860.0 |  1 |  297 |  1 | 0.00%
|  p1.0n40s1 |  23.994 |  23.789 |  3,019 |  2,313 |  335.4 |  9 |  534 |  3 | 0.85%
|  p1.0n40s2 |  22.401 |  22.401 |  2,376 |  2,376 |  2,376.0 |  1 |  262 |  1 | 0.00%
|  p1.0n40s3 |  21.942 |  21.942 |  3,300 |  3,300 |  3,300.0 |  1 |  274 |  1 | 0.00%
|  p1.0n40s4 |  18.342 |  18.228 |  3,545 |  1,655 |  154.1 |  23 |  854 |  4 | 0.62%
Correlation valuation |  p0.6pS0.6n35s0 |  224.0 |  224 |  256 |  256 |  256.0 |  1 |  433 |  1 | 0.00%
|  p0.6pS0.6n35s1 |  243.0 |  243 |  458 |  350 |  50.8 |  9 |  688 |  2 | 0.00%
|  p0.6pS0.6n35s2 |  239.0 |  239 |  451 |  451 |  451.0 |  1 |  671 |  1 | 0.00%
|  p0.6pS0.6n35s3 |  241.8 |  241 |  711 |  222 |  28.4 |  25 |  1,684 |  1 | 0.33%
|  p0.6pS0.6n35s4 |  248.5 |  248 |  465 |  392 |  155 |  3 |  705 |  1 | 0.20%
|  p0.6pS0.6n40s0 |  292.0 |  292 |  811 |  811 |  811.0 |  1 |  538 |  1 | 0.00%
|  p0.6pS0.6n40s1 |  298.0 |  298 |  1,778 |  1,778 |  1,778.0 |  1 |  1,130 |  1 | 0.00%
|  p0.6pS0.6n40s2 |  289.0 |  289 |  1,735 |  1,735 |  1,735.0 |  1 |  596 |  1 | 0.00%
|  p0.6pS0.6n40s3 |  311.0 |  311 |  938 |  938 |  938.0 |  1 |  500 |  1 | 0.00%
|  p0.6pS0.6n40s4 |  313.5 |  313 |  4,542 |  3,062 |  1,514 |  3 |  1,565 |  1 | 0.15%
|  p0.8pS0.6n35s0 |  304.0 |  304 |  1,714 |  1,277 |  342.8 |  5 |  829 |  3 | 0.00%
|  p0.8pS0.6n35s1 |  294.5 |  294 |  1,746 |  834 |  116.4 |  15 |  1,614 |  2 | 0.16%
|  p0.8pS0.6n35s2 |  299.0 |  299 |  1,670 |  1,670 |  1,670.0 |  1 |  704 |  1 | 0.00%
|  p0.8pS0.6n35s3 |  307.0 |  307 |  1,070 |  1,070 |  1,070.0 |  1 |  567 |  1 | 0.00%
|  p0.8pS0.6n35s4 |  298.0 |  298 |  1,137 |  1,137 |  1,137.0 |  1 |  422 |  1 | 0.00%
|  p0.8pS0.6n40s0 |  395.0 |  395 |  4,271 |  4,271 |  4,271.0 |  1 |  839 |  1 | 0.00%
|  p0.8pS0.6n40s1 |  380.5 |  380 |  6,643 |  4,522 |  2,214.0 |  3 |  1,500 |  1 | 0.13%
|  p0.8pS0.6n40s2 |  385.0 |  385 |  6,402 |  6,402 |  6,402.0 |  1 |  850 |  1 | 0.00%
|  p0.8pS0.6n40s3 |  394.0 |  394 |  9,791 |  9,791 |  9,791.0 |  1 |  897 |  1 | 0.00%
|  p0.8pS0.6n40s4 |  387.0 |  387 |  7,278 |  7,278 |  7,278 |  1 |  734 |  1 | 0.00%

We slightly generalized the use of the correlation function by supposing that
some edges don’t have signs at all, that is, allowing an incomplete underlying
graph. The random graphs were generated using the same model and an edge has a
”$+$” or a ”$-$” sign with a prescribed probability. The name of a benchmark
file for this valuation indicates the probability that an edge belongs to the
graph (”p”), the probability of the plus sign (”pSign”), the number of agents
(”n”), and finally the number of the instance (”s”).

Table 1 shows the results of our numerical tests; it contains the optimum
value in the root (LP problem), the (most of the time) optimum value of the
ILP problem found by branch and bound algorithm, the overall running time, the
running time for finding and solving the LP problem in the root, the running
time per node, the number of nodes of the branching tree, the total number of
variables added during the execution, the number of nodes fathomed by
integrality, and the gap between the LP and the ILP optimum values.

Our algorithm is one of the first to present competitive results for more than
$20$ agents for quite general valuation functions, using only a regular PC.
For half of the tests the root LP problem has an integer solution and for all
of them the gap is less than one percent. The branching trees (with few
exceptions) are small and most of the running time is concentrated in the root
and in the nodes corresponding to variables set to zero or variables set to
one but associated with small sized coalitions. Conducting the branching rule
to find medium or large sized coalitions variables seems to be of secondary
interest since the resulting trees are small in size.

We listed only the results using the quadratic models for solving the sub-
problem, however the linear models have a good potential which deserves future
exploring work (Lagrangian relaxation, special treatment of the high
degeneracy etc). Future work should also be directed towards getting rid of
the main burden of our algorithm: find specific methods for solving QKPf.
Another direction of future research could be a cut generation method (finding
inequalities defining facets of a subjacent polyhedral relaxation) that
tightens the LP relaxations towards a branch and cut algorithm.

##  Appendix A Appendix: The coordination valuation

Let $G=(V,E)$ be graph and $C\subseteq V$ be a coalition, we define

$\displaystyle\displaystyle n_{i}(C)=|\\{jk\in E\>:\>j\in C,k\notin C,ij,ik\in
E\\}|.$

The coordination coalition valuation function ([13]) is

$v:2^{V}\to{\mathbb{R}},v(C)=\sum_{i\in C}n_{i}(C),\forall C\subseteq V.$

This function accounts for all cliques of three agents, two of them being
inside the coalition, while the third is outside. It is designed to offer
coalitions that include agents that have common neighbours from outside. We
will generalize this definition by including the weight of the $3$-clique as
the sum of the weights on its edges. Hence, we consider first
$w:E\to{\mathbb{R}}$ to be an weight on the edges of $G$. Second, if
${\mathbf{v}}$ is the characteristic vector for the generic coalition $C$,
then

$\displaystyle\begin{array}[]{l}\displaystyle
v(C)=\sum_{i=1}^{n}v_{i}\sum_{j=1}^{n}w_{ij}v_{j}\sum_{k=1}^{n}w_{ik}w_{jk}(1-v_{k}).\end{array}$
(99)

We get the coordination valuation function by taking $w$ to be the edge
characteristic function, i. e., $w_{ij}=1$, if $ij\in E$, and $0$ otherwise.
The sub-problem (33) becomes:

$\displaystyle\begin{array}[]{rcl}&max&\displaystyle\left(\sum_{i=1}^{n}\sum_{j=1}^{n}\sum_{k=1}^{n}w_{ij}w_{ik}w_{jk}v_{i}v_{j}(1-v_{k})+\sum_{i\in
V}\pi_{i}v_{i}\right)\end{array}$ (101)
$\displaystyle\begin{array}[]{rcl}\displaystyle\sum_{i\in
C_{j}}v_{i}+\sum_{i\notin C_{j}}(1-v_{i})&\geq&1,\forall j\in J\end{array}$
(103) $\displaystyle\begin{array}[]{rcl}\displaystyle
v_{i}&\in&\\{0,1\\},\forall i\in V\end{array}$ (105)

###  A.1 Reducing the degree of the objective function

Using the same method of Glover from ([6]), we will transfom the cubic
objective function from (101) to a quadratic one by adding only $n$ new
unknowns and $4n$ new (quadratic) constraints. We introduce first the
variables

$\displaystyle\begin{array}[]{l}\displaystyle
z_{i}=v_{i}\sum_{j=1}^{n}v_{j}\left(\sum_{k=1}^{n}w_{ij}w_{ik}w_{jk}(1-v_{k})\right),\forall
i\in V.\end{array}$ (107)

Second, add the constraints

$\displaystyle\begin{array}[]{l}Y_{i}^{-}v_{i}\leq z_{i}\leq
Y_{i}^{+}v_{i},\forall i\in V\\\
\displaystyle\sum_{j=1}^{n}v_{j}\left(\sum_{k=1}^{n}w_{ij}w_{ik}w_{jk}(1-v_{k})\right)-Y_{i}^{+}(1-v_{i})\leq
z_{i},\forall i\in V\\\
\displaystyle\sum_{j=1}^{n}v_{j}\left(\sum_{k=1}^{n}w_{ij}w_{ik}w_{jk}(1-v_{k})\right)-Y_{i}^{-}(1-v_{i})\geq
z_{i},\forall i\in V,\end{array}$ (111)

where $Y_{i}^{-}$ and $Y_{i}^{+}$ are a lower and, respectively, an upper
bound for the expression

$\displaystyle\begin{array}[]{l}\displaystyle\sum_{j=1}^{n}v_{j}\left(\sum_{k=1}^{n}w_{ij}w_{ik}w_{jk}(1-v_{k})\right).\end{array}$
(113)

$Y_{i}^{-}$ and $Y_{i}^{+}$ could be

$\displaystyle\begin{array}[]{l}Y_{i}^{-}=\displaystyle\left(\sum_{j,k:w_{ij}w_{ik}w_{jk}<0}w_{ij}w_{ik}w_{jk}\right)\left/2\right.\\\
Y_{i}^{+}=\displaystyle\left(\sum_{j,k:w_{ij}w_{ik}w_{jk}>0}w_{ij}w_{ik}w_{jk}\right)\left/2\right.\end{array}$
(116)

The problem (101) - (105) becomes

$\displaystyle\begin{array}[]{rcl}&max&\displaystyle\left(\sum_{i=1}^{n}\pi_{i}v_{i}+\sum_{i=1}^{n}z_{i}\right)\end{array}$
(118) $\displaystyle\begin{array}[]{rcl}\displaystyle\sum_{i\in
C_{j}}v_{i}+\sum_{i\notin C_{j}}(1-v_{i})&\geq&1,\forall j\in J,\end{array}$
(120)
$\displaystyle\begin{array}[]{rcl}\displaystyle\sum_{j=1}^{n}v_{j}\left(\sum_{k=1}^{n}w_{ij}w_{ik}w_{jk}(1-v_{k})\right)+Y_{i}^{+}v_{i}-z_{i}\leq
Y_{i}^{+},\forall i\in V,\end{array}$ (122)
$\displaystyle\begin{array}[]{rcl}\displaystyle\sum_{j=1}^{n}v_{j}\left(\sum_{k=1}^{n}w_{ij}w_{ik}w_{jk}(1-v_{k})\right)+Y_{i}^{-}v_{i}-z_{i}\geq
Y_{i}^{-},\forall i\in V,\end{array}$ (124)
$\displaystyle\begin{array}[]{rcl}Y_{i}^{-}v_{i}-z_{i}\leq 0,\forall i\in
V,\end{array}$ (126)
$\displaystyle\begin{array}[]{rcl}Y_{i}^{+}v_{i}-z_{i}&\geq&0,\forall i\in
V,\end{array}$ (128) $\displaystyle\begin{array}[]{rcl}\displaystyle
v_{i}&\in&\\{0,1\\},\forall i\in V,\end{array}$ (130)
$\displaystyle\begin{array}[]{rcl}\displaystyle z_{i}&\in&{\mathbb{R}},\forall
i\in V.\end{array}$ (132)

(For the coordination valuation function we can simply take $Y_{i}^{-}=0$ and
$Y_{i}^{+}=n^{2}$, for all $i\in V$.)

This last problem is a Mixed Integer Quadratically Constrained Programming
(MIQCP) problem that could also be solved using a mathematical optimization
solver. Another way of addressing this problem is to apply again the method of
Glover to the expressions

$\displaystyle\begin{array}[]{l}z_{ij}=\displaystyle
v_{j}\left(\sum_{k=1}^{n}w_{ij}w_{ik}w_{jk}(1-v_{k})\right),\end{array}$ (134)

which will add $n^{2}$ new variables and $4n^{2}$ new constraints to the
problem.

###  A.2 Branch and Bound version for the coordination valuation function

Suppose as above that the current node in the branching tree has a subset of
already covered agents, $\displaystyle U=\displaystyle\bigcup_{j\in
J:x_{j}=1}C_{j}$, and $U^{\prime}=V\setminus U$. For this valuation function
the sub-problem (33) becomes:

$\displaystyle\begin{array}[]{rcl}&max&\displaystyle\left[\sum_{i\in
U^{\prime}}\sum_{j\in U^{\prime}}\left(\sum_{k\in
U^{\prime}}w_{ij}w_{ik}w_{jk}v_{i}v_{j}(1-v_{k})+\sum_{k\in
U}w_{ij}w_{ik}w_{jk}v_{i}v_{j}\right)+\right.\\\
&&\displaystyle\left.+\sum_{i\in U^{\prime}}\pi_{i}v_{i}\right]\end{array}$
(137) $\displaystyle\begin{array}[]{rcl}\displaystyle\sum_{i\in
C_{j}}v_{i}+\sum_{i\notin C_{j}}(1-v_{i})&\geq&1,\forall j\in J,\end{array}$
(139) $\displaystyle\begin{array}[]{rcl}\displaystyle
v_{i}&\in&\\{0,1\\},\forall i\in U^{\prime}.\end{array}$ (141)

Using the above degree reduction procedure for the objective function we
introduce the variables

$\displaystyle\begin{array}[]{l}\displaystyle z_{i}=v_{i}\sum_{j\in
U^{\prime}}v_{j}\left(\sum_{k\in
U^{\prime}}w_{ij}w_{ik}w_{jk}(1-v_{k})+\sum_{k\in
U}w_{ij}w_{ik}w_{jk}\right),\forall i\in U^{\prime}.\end{array}$ (143)

Second, add the constraints

$\displaystyle\begin{array}[]{l}Y_{i}^{-}v_{i}\leq z_{i}\leq
Y_{i}^{+}v_{i},\forall i\in V\\\ \displaystyle\sum_{j\in
U^{\prime}}v_{j}\left(\sum_{k\in
U^{\prime}}w_{ij}w_{ik}w_{jk}(1-v_{k})+\sum_{k\in
U}w_{ij}w_{ik}w_{jk}\right)-Y_{i}^{+}(1-v_{i})\leq z_{i},\forall i\in
U^{\prime}\\\ \displaystyle\sum_{j\in U^{\prime}}v_{j}\left(\sum_{k\in
U^{\prime}}w_{ij}w_{ik}w_{jk}(1-v_{k})+\sum_{k\in
U}w_{ij}w_{ik}w_{jk}\right)-Y_{i}^{-}(1-v_{i})\geq z_{i},\forall i\in
U^{\prime},\end{array}$ (147)

where $Y_{i}^{-}$ and $Y_{i}^{+}$ are a lower and, respectively, an upper
bound for the expression

$\displaystyle\begin{array}[]{l}\displaystyle\sum_{j\in
U^{\prime}}v_{j}\left(\sum_{k\in
U^{\prime}}w_{ij}w_{ik}w_{jk}(1-v_{k})+\sum_{k\in
U}w_{ij}w_{ik}w_{jk}\right).\end{array}$ (149)

$Y_{i}^{-}$ and $Y_{i}^{+}$ could be

$\displaystyle\begin{array}[]{l}Y_{i}^{-}=\displaystyle\left(\sum_{j,k\in
U^{\prime}:w_{ij}w_{ik}w_{jk}<0}w_{ij}w_{ik}w_{jk}\right)\left/2\right.+\sum_{j\in
U^{\prime},k\in U:w_{ij}w_{ik}w_{jk}<0}w_{ij}w_{ik}w_{jk}\\\
Y_{i}^{+}=\displaystyle\left(\sum_{j,k\in
U^{\prime}:w_{ij}w_{ik}w_{jk}>0}w_{ij}w_{ik}w_{jk}\right)\left/2\right.+\sum_{j\in
U^{\prime},k\in U:w_{ij}w_{ik}w_{jk}>0}w_{ij}w_{ik}w_{jk}\end{array}$ (152)

The problem (101) - (105) becomes

$\displaystyle\begin{array}[]{rcl}&max&\displaystyle\left(\sum_{i\in
U^{\prime}}z_{i}+\sum_{i\in W^{\prime}}\pi_{i}v_{i}\right)\end{array}$ (154)
$\displaystyle\begin{array}[]{rcl}\displaystyle\sum_{i\in
C_{j}}v_{i}+\sum_{i\notin C_{j}}(1-v_{i})&\geq&1,\forall j\in J,\end{array}$
(156) $\displaystyle\begin{array}[]{rcl}\displaystyle\sum_{j\in
U^{\prime}}v_{j}\left(\sum_{k\in
U^{\prime}}w_{ij}w_{ik}w_{jk}(1-v_{k})+\sum_{k\in
U}w_{ij}w_{ik}w_{jk}\right)+Y_{i}^{+}v_{i}-z_{i}\leq Y_{i}^{+},\forall i\in
U^{\prime},\end{array}$ (158)
$\displaystyle\begin{array}[]{rcl}\displaystyle\sum_{j\in
U^{\prime}}v_{j}\left(\sum_{k\in
U^{\prime}}w_{ij}w_{ik}w_{jk}(1-v_{k})+\sum_{k\in
U}w_{ij}w_{ik}w_{jk}\right)+Y_{i}^{-}v_{i}-z_{i}\geq Y_{i}^{-},\forall i\in
U^{\prime},\end{array}$ (160)
$\displaystyle\begin{array}[]{rcl}Y_{i}^{-}v_{i}-z_{i}\leq 0,\forall i\in
U^{\prime},\end{array}$ (162)
$\displaystyle\begin{array}[]{rcl}Y_{i}^{+}v_{i}-z_{i}&\geq&0,\forall i\in
U^{\prime},\end{array}$ (164) $\displaystyle\begin{array}[]{rcl}\displaystyle
v_{i}&\in&\\{0,1\\},\forall i\in U^{\prime},\end{array}$ (166)
$\displaystyle\begin{array}[]{rcl}\displaystyle z_{i}&\in&{\mathbb{R}},\forall
i\in U^{\prime}.\end{array}$ (168)

(For the coordination valuation function we can simply take $Y_{i}^{-}=0$ and
$Y_{i}^{+}=|U^{\prime}|^{2}/2+|U^{\prime}|\cdot|U|\leq 3n^{2}/4$, for all
$i\in U^{\prime}$.)

## References

  * [1] Bansal, N., Blum, A., Chawla, S.,  Correlation Clustering, Machine Learning, 56, 89–113, 2004.
  * [2] Bistaffa, F., Farinelli, A., Cerquides, J., Rodriguez-Aguilar, J., Ramchurn, S. D.,  Anytime coalition structure generation on synergy graphs, Proceedings of AAMAS’14,13–20, 2014.
  * [3] Deng, X., Papadimitriou., C.,  On the complexity of cooperative solution concepts, Mathematics of Operations Research, 19 (2), 257–266, 1994.
  * [4] Flammini, M., Monaco, G., Moscardelli, L., Shalom, M., Zaks, S.,  Online Coalition Structure Generation in Graph Games,Proceedings of AAMAS’18, 1353–1361, 2018.
  * [5] Gallo, G., Hammer, P.L., Simeone, B.,  Quadratic Knapsack Problems, Mathematical Programming, 12 (2), 132–149, 1980.
  * [6] Glover, F, Wolsey, E.,  Improved linear integer programming formulations of nonlinear integer problems, Amanagement Science, 22 (4), 180–182, 1975.
  * [7] Hoffman, K. L.., Padberg, M.,  Solving Airline Crew Scheduling Problems by Branch-and-Cut, Management Science, 6 (39), 657–682, 1993.
  * [8] Pisinger, D.,  The Quadratic Knapsack Problem - a Survey, Discrete Applied Mathematics, 155, 623–648, 2007.
  * [9] Rahwan, T., Ramchurn, S. D., Dang, V. D., Giovanucci, A., Jennings, N. R.,  Anytime optimal coalition structure generation, Proceedings of AAAI’07, 1184–1190, 2007.
  * [10] Rahwan, T., Michalak, T. P., Wooldridge, M., Jennings, N. R.,  Coalition structure generation: A survey, Artificial Intelligence, 229, 139–174, 2015.
  * [11] Sandholm, T., Larson, K., Andersson, M., Shehory, O., Tohme, F.,  Coalition structure generation with worst case guarantees, Artificial Intelligence, 1–2, 209–238, 2009.
  * [12] Ueda, S., Iwasaki, A., Conitzer, V., Ohta, N., Sakurai, Y., Yakoo, M.,  Coalition structure generation in cooperative games with compact representations, Autonomous Agents and Multi-Agent Systems, 32, 503–533, 2018.
  * [13] Voice, T., Polukarov, M., Jennings, N. R.,  Coalition Structure Generation over Graphs, Journal of Artificial Intelligence Research, 45, 165–196, 2012.
  * [14] Voice, T., Ramchurn, S., Jennings, N. R.,  On coalition formation with sparse synergies, Proceedings of AAMAS’12, 45, 223–230, 2012.


# Bridging Walking and Slithering – Stokesian Locomotion

Shai Revzen, Dan Zhao, Glenna Clifton, Brian Bittner, Nick Gravish

Presented at Dynamic Walking on May 20th 2021

##  I Introduction

Both legged locomotion and slithering motions typically utilize periodic gaits
– repeating cycles of body shape change that produce a net motion through the
world. Legged locomotion can be viewed from the perspective of piece-wise
contact constraint formation and removal. Slithering and low Reynolds number
swimming operate under continuous constraints of force balance, wherein
dissipation removes the ability to accumulate momentum. Here we discuss how to
bridge the gap between these domains of motion, thereby, among other benefits,
producing models for the space of legged locomotion _with slipping_. The
connective fabric is the use of a “Stokesian”, or “local connection” model of
the form $g^{-1}\dot{g}=A(r)\dot{r}$ where $g\in\mathsf{SE}(3)$ is the body
frame location and $r\in\mathcal{Q}$ is the body shape.

##  II Holonomic Motion

Leading up to locomotion, examine the ability of a robotic manipulator to move
and orient its own center of mass111We tacitly assume a CoM attached “body
frame”; the proper selection of such a frame is non-trivial [2].. There is
typically a map $F:\mathcal{Q}\to\mathsf{SE}(3)$ from the configuration
parameters describing the shape of the robot $r\in\mathcal{Q}$, to the body
frame in $g=F(r)\in\mathsf{SE}(3)$. Here the constraints are “holonomic”, i.e.
they arise from the “whole” function $F$, and as a consequence no periodic
motion in $r$ can lead to a net motion through the world. When $r$ returns to
its initial value, the $g$ also returns. The constraint in this case can be
studied using the Jacobian $J_{F}$ of $F$. We can construct a body velocity
map as a “local connection” $A_{F}(r):=(F(r))^{-1}J_{F}(r)$ and thus
$A_{F}(r)\dot{r}\in\mathfrak{se}(3)$ is a body velocity222We assumed notation
of a matrix Lie group representation for $\mathsf{SE}(3)$. .

##  III Piece-wise holonomic motion

Typically, legged locomotion is modeled assuming non-slip contacts with the
environment. For simplicity, we restrict our treatment to sets of contact
constraints that fully determine the body frame location relative to the feet,
i.e. are not under-actuated. These could be three or more non-colinear pin-
joint foot contacts, two foot contacts at least one of which can bear a moment
around the line connecting them, or one foot which can bear all three moments.

It is important to note that even for a simple case where which feet are in
contact is uniquely determined by the body shape, $g$ can no longer be written
as a function of $r$. For every continuous set of contacts $c\in C$ out of a
finite set of contact combinations $C$, there exists an
$F[c]:\mathcal{Q}\to\mathsf{SE}(3)$. If any one of the contact combinations is
maintained no motion is possible for the same reason above. However motion is
possible, because the change in contacts involves a jump discontinuity in $g$
– the $F[c]$ do not meet continuously at the switching of $c$, i.e. in general
$F[c_{1}](r)\neq F[c_{2}](r)$ even for $r$ values where $c_{1}$ can switch to
$c_{2}$.

Here the local connection allows for a concise and precise representation
using $A_{F[c]}(r)$. Like $F[c]$ this function has a jump discontinuity when
$c$ changes, but it allow us to model motion using the ODE
$g^{-1}\dot{g}=A_{F[c]}(r)\dot{r}$, which produces continuous solutions in $g$
– the body frame.

##  IV Slithering Motion

For slithering and low Reynolds number (so called “Stokes Fluid”) swimming, it
is well established that a similar mathematical formalism
$g^{-1}\dot{g}=A(r)\dot{r}$ applies. However here $A(r)$ is a fundamental
object, and nowhere is it derived from an underlying Jacobian. This expresses
itself in the fact that most Stokesian systems will move the body frame with
almost every cyclic shape change, unlike piecewise holonomic systems which
require a cycle that spans multiple of the individually holonomic pieces.
$A(r)$ could be approximated piece-wise with pieces that come from holonomic
systems; in this sense slithering could arise as a limit of walking with an
infinite number of legs.

##  V Walking with slipping

Consider the removal of the non-slip requirements that typically produce
holonomic motion. With slipping $A(r)$ for a fixed set of contacts $c$ is no
longer derived from a Jacobian of some underlying $F[c]$. Instead, it is
derived from a set of implicit equations tying together body frame velocity
and shape change speed $\dot{r}$, i.e. $G(r;g^{-1}\dot{g},\dot{r})=0$ If the
$G$ is linear in $\dot{r}$ and in $g^{-1}\dot{g}$, the resulting system also
takes the form $g^{-1}\dot{g}=A(r)\dot{r}$.

##  VI Summary

The common mathematical formalism $g^{-1}\dot{g}=A(r)\dot{r}$ extends across
all four types of motion – holonomic, piece-wise holonomic, slithering, and
walking with slipping. It is the special case of the “reconstruction equation”
of [1] for the case where momentum can be ignored (see also [3]).

## References

  * [1] Bloch, Marsden and Zenkov  “Nonholonomic dynamics”  In _Not. AMS_ 52.3, 2005, pp. 320–329
  * [2] Hatton and Choset  “Geometric motion planning: The local connection, Stokes theorem, and the importance of coordinate choice”  In _IJRR_ 30.8 SAGE Publications Sage UK: London, England, 2011, pp. 988–1014
  * [3] Kvalheim, Bittner and Revzen  “Gait modeling and optimization for the perturbed Stokes regime”  In _Nonlinear Dynamics_ , 2019  DOI: 10.1007/s11071-019-05121-3