Showing posts with label Optimisation. Show all posts
Showing posts with label Optimisation. Show all posts

Sunday, 11 May 2014

Wizards, Giants, Linear Programming and Sage

I've been playing a game on my phone for a couple of months now called Clash of clans. It's one of those freemium games that tells me to do something on my phone every now and then:




In essence during the game you build your base with buildings to obtain resources, create troops and defend against attacks.

Here's a screen shot of my base at the moment:


As well as building your base, you can also put together a raiding party and attacking other bases.
In this post, I'll describe the use of a mathematical technique called linear programming so as to get the best combination of Giants, Wizards etc in a raiding party.

To train a warrior costs Elixir (a resource that you mine and/or plunder) but also costs space in your army camps.

Here are the various warriors available to me at the moment (as you build better buildings you get more warriors etc...).

1. Barbarian (Damage per second: 18, Hitpoints: 78, Housing Space: 1, Elixir Cost: 80):


2. Archer (Damage per second: 16, Hitpoints: 33, Housing Space: 1, Elixir Cost: 160)


3. Goblin (Damage per second: 19, Hitpoints: 36, Housing Space: 1, Elixir Cost: 60)


4. Giant (Damage per second: 24, Hitpoints: 520, Housing Space: 5, Elixir Cost: 2000) 


 5. Wall Breaker (Damage per second: 32, Hitpoints: 35, Housing Space: 2, Elixir Cost: 2500)


6. Balloon (Damage per second: 72, Hitpoints: 280, Housing Space: 5, Elixir Cost: 3500)


7. Wizards (Damage per second: 125, Hitpoints: 130, Housing Space: 4, Elixir Cost: 3000)

8. Healer (Damage per second: NA, Hitpoints: 600, Housing Space: 14, Elixir Cost: 6000)


9. Dragon (Damage per second: 140, Hitpoints: 1900, Housing Space: 20, Elixir Cost: 25000)




To summarise the information about the players I'm going to put the costs in a vector \(c\) and the housing space in a vector \(h\) (indexed using the same order as the images above):

\[c=(80,160,60,2000,2500,3500,3000,6000,25000)\]
\[h=(1,1,1,5,2,5,4,14,20)\]

Currently I have a capacity of 50 housing spaces per camp and I have 4 camps so I have an upper bound on the total amount of housing space of \(K=200\).

If I let the vector \(x\) denote a raiding party so that \(x_i\) is the number of warriors of type \(\i) (using the ordering above) then we have the corresponding housing space constraint:

\[\sum_{i=1}^9h_ix_i\leq K\]

I can now try and do two things:
  1. Get a full set of camps by minimising the overall cost.
  2. Get the most powerful raiding party that can fit in my camps.
I will solve both of these problems using what is called Integer Linear Programming. The integer part of that is due to the fact that I must have whole parts of warriors.

  1. Minimising my costs

The cost of a given rading party \(x\) is given by:

\[C(x) = \sum_{i=1}^9c_ix_i \]

Thus to minimise this cost and fill the camps we need to solve the following minimisation problem:

\[\text{minimise }C(x) = \sum_{i=1}^9c_ix_i\]

such that:

\[H(x)=\sum_{i=1}^9h_ix_i= K\]
and
\[x_i\in\mathbb{Z}_{\geq0}\]

How do we solve this? I've blogged about a similar type of optimisation problem that I use to schedule my class assessment and for that I used the awesome sagemath. You can read that blog post here. Here I will do the same thing.

The following Sage code is what is needed to create the linear program and also obtain the solution:

# Set up parameters:
K = 200  # Set the upper bound
c = [80,160,60,2000,2500,3500,3000,6000,25000] # The the cost vector
h = [1,1,1,5,2,5,4,14,20]  # Set the 'housing space' constraints vector
n = len(c)  # Get the number of warriors

C = lambda x: sum([x[i]*c[i] for i in range(n)])  # Define the cost function
H = lambda x: sum([x[i]*h[i] for i in range(n)])  # Define the housing capacity function

# Create and solve the optimisation problem
p = MixedIntegerLinearProgram(maximization=False)  # Create an optimisation problem with minimization as the goal
x = p.new_variable(integer=True)  # Create a variable with the requirement that it must be integer
p.add_constraint(H(x)==K)  # Set the housing constranit
p.set_objective(C(x))  # Set the objective function
p.show() # Show the optimisation problem
p.solve()  # Solve the optimisation problem

# Solve the optimisation problem

print 'Has solution:' 
for i, v in p.get_values(x).iteritems(): 
    print 'x_%s = %s' % (i, int(round(v))) 
print 'Housing spaces used: %s' % H(p.get_values(x))

The output is (disappointing but expected): \(x=(0,0,200,0,0,0,0,0)\) which implies that I should have 200 Goblins.

If you're familiar with the game this is not a very sensible way to go as Goblins are really only good at getting loot from defenceless bases.

We can add a couple more constraints so that we don't have too many goblins (let's limit it to 10) and also so that we have at least one healer:

p = MixedIntegerLinearProgram(maximization=False)  # Create an optimisation problem with minimization as the goal
x = p.new_variable(integer=True)  # Create a variable with the requirement that it must be integer
p.add_constraint(H(x)==K)
p.add_constraint(x[2]<=10)  # Don't want more than 10 goblins
p.add_constraint(x[7]>=1)  # Want 1 healer at least
p.set_objective(C(x))
p.show()
p.solve()
print 'Has solution:'
for i, v in p.get_values(x).iteritems():
    print 'x_%s = %s' % (i, int(round(v)))
print 'Housing spaces used: %s' % H(p.get_values(x))

This gives us \(x=(176,0,10,0,0,0,0,1,0)\).

This is again not very satisfying as we're basically filling our raiding party with the cheapest options. That was the purpose of this approach though. 

The next thing we'll look at is how to actually get a good raiding part.

  1. Getting the 'best' raiding party.

As you can see on the images, every warrior is different: some fly, some are better against buildings, some have more hitpoints than others etc...

One way of getting the 'best' raiding party could be to maximise the damage per second of the raiding party. To do this, let's define the following vector (which is just the damage per second that can be read off of the images above):

\[d=(18, 16, 19, 24, 32, 72, 125, 0, 140)\]

Now our maximisation problem can be written as:

\[\text{maximise } \sum_{i=1}^9d_ix_i\]

such that:

\[H(x)=\sum_{i=1}^9h_ix_i\leq K\]
and
\[x_i\in\mathbb{Z}_{\geq0}\]

Here's the Sage code to do this:

damagepersecond = [18, 16, 19, 24, 32, 72, 125, 0, 140]
C = lambda x: sum([x[i]*damagepersecond[i] for i in range(n)])

p = MixedIntegerLinearProgram()  # Create an optimisation problem with minimization as the goal
x = p.new_variable(integer=True)  # Create a variable with the requirement that it must be integer
p.add_constraint(H(x)<=K)
p.add_constraint(x[2]<=10)  # Don't want more than 10 goblins
p.add_constraint(x[7]>=1)  # Want 1 healer at least
p.set_objective(C(x))
p.show()
p.solve()
print 'Has solution:'
for i, v in p.get_values(x).iteritems():
    print 'x_%s = %s' % (i, int(round(v)))
    print 'Housing spaces used: %s' % H(p.get_values(x))

The result of the above is \(x=(0,0,2,0,0,0,46,1,0)\). I.e take 2 goblins, 46 wizards and 1 healer.
The problem with this raiding party is that wizards are actually quite fragile (archers towers and cannons can pretty much just pick them off).

So perhaps instead of maximising the damage per second we could maximise the total hit points of the party.

To do this we need to define the vector of hitpoints:

\[p=(78, 33, 36, 520, 35, 280, 130, 600, 1900)\]

Now our maximisation problem can be written as:

\[\text{maximise } \sum_{i=1}^9p_ix_i\]

such that:

\[H(x)=\sum_{i=1}^9h_ix_i\leq K\]
and
\[x_i\in\mathbb{Z}_{\geq0}\]

This gives us \(x=(0,0,0,0,93,0,0,1,0)\): a raiding part of wall breakers.

The problem with this is that wall breakers are good for just that: breaking walls. They carry a bomb up to a wall and then they explode. So let's add some more constraints to our optimisation problem so as to get something a bit more useful.
  • Let's make sure that we always have at least 20 housing capacity worth of distance attackers (Archers or Wizards):
\[x_1h_1+x_6h_6 \geq 20 \]
  • Let's ensure that we have at least 1 air attacker (so that once the air defenses are out of the way we can finish off the base):
\[x_5+x_8\geq1\]
  • Let's ensure that we always have some wall breakers to clear the walls for our ground units (1 per 5 giants and 1 per 20 barbarians say):
\[5x_4\geq x_3\]

\[20x_4\geq x_0\]
The code to solve this is given here:

p = MixedIntegerLinearProgram()  # Create an optimisation problem with minimization as the goal
x = p.new_variable(integer=True)  # Create a variable with the requirement that it must be integer
p.add_constraint(H(x) <= K)
p.add_constraint(x[2] <= 10)  # Don't want more than 10 goblins
p.add_constraint(x[7] >= 1)  # Want 1 healer at least
p.add_constraint(x[1] * h[1] + x[6] * h[6]>=20)  # Want at least 20 capacity worth of distance attackers
p.add_constraint(x[5] + x[8]>=1)  # Want at least 1 of air attackers
p.add_constraint(5 * x[4] >= x[3])  # Want at least 1 wall breaker per 5 giants
p.add_constraint(20 * x[4] >= x[0])  # Want at least 1 wall breaker per 20 barbarians
p.set_objective(C(x))
p.show()
p.solve()
print 'Has solution:'
for i, v in p.get_values(x).iteritems():
    print 'x_%s = %s' % (i, int(round(v)))
print 'Housing spaces used: %s' % H(p.get_values(x))

This gives \(x=(1,20,0,23,5,0,0,1,2)\): so take 1 Barbarian, 20 Archers, 23 Giants, 5 Wall Breakers, 1 Healer and 2 Dragons.

The above would give us the most hit points but now we're missing out on the damage power of wizards.

Perhaps we could try and balance both hitpoints damage.

In essence we have what's called a multi objective optimisation problem. One approach to solve this is to simply weight both objectives.

We're going to use \(\alpha\) to denote the importance of the hitpoints and \(1-\alpha\) the importance of the damage.

Our optimisation function then becomes:

\[\alpha\sum_{i=1}^9p_ix_i+(1-\alpha)\sum_{i=1}^9d_ix_i\]

Using the same constraints as before and taking \(\alpha=.5\) (so that both damage and hit points have equal weighting) we get \(x=(11,0,0,25,5,0,5,1,1)\). Thus we should be taking in to battle: 11 Barbarians, 25 Giants, 5 Wall Breakers, 5 Wizards, 1 Healer and 1 Dragon.

This feels much more balanced than anything we've had so far and not far off what I've been using as my own raiding party.

Given that everything is nicely coded we can quickly obtain the make up of the optimal raiding party for a range of values of \(\alpha\):



We see that as \(\alpha\) increases (so that we care more about hit points) the number of Wizards we should use slowly decreases (at one point the constraint for ranged attacking units kicks in and we replace the Wizards with Archers).

You could easily tweak the constraints to perhaps ensure your party always has some Goblins or otherwise. I think this is a nice example where a mathematical model can be used to solve a problem but at the same time must not be taken far from the problem itself. There might be some players who only choose to attack bases with a particular type of defence so they can tweak/add constraints to the problem accordingly.

Feel free to play around with the code in the Sage cell below:


The mathematical techniques used in this post rely on Linear Algebra, Higher Dimensional Geometry and obviously the ability to write code :)

(On another note I'm still not in a clan on the game so if anyone has any space in a clan let me know!)

Thursday, 27 March 2014

Scheduling group presentations using Graph Theory and Sage.

Yesterday (2014-03-26) I spoke at the Embedded Enterprise Exchange about my new first year module. I used a Hangout on air during the talk and you can see it here:

That's not what I'm going to talk about here though. The course I talked baout ends with all the 'companies' (groups of 4 students) giving a $\approx 25$ minute talk.

I need to schedule 40ish talks and this needs to fit around the student availability as well as my own. In this post I'll describe how I did this using a combination of Doodle, +Sage Mathematical Software System and Graph Theory.

The beginning of this post will be some of the boring details but towards the end I start talking about the mathematics (so feel free to skip to there...).

First of all I needed to know my students availability.

For this I simply used Doodle: https://doodle.com. I kind of use Doodle for every meeting I have to schedule (they also offer a cool tool that lets you show your availability so students/colleagues have an indication of when I might be able to meet with them).

Here's a screenshot of the responses:


You can't really see anything there as I had to zoom out a lot to grab the whole picture. Doodle allows you to download the information for any given poll in .xls format so I could relatively easily obtain the biadjacency matrix $M$ for my problem. Where $M_{ij}$ is 1 if group $i$ is available for schedule slot $j$ and 0 otherwise.


The mathematics and code needed.

Once I've got a .csv file (by tweaking the .xls file) of the biadjacency matrix I import that in to +Sage Mathematical Software System  and convert it to an instance of the `Matrix` class using the following:

import csv 
f=open(DATA+'availabilitymatrix','r') 
data = [[int(j) for j in row] for row in csv.reader(f)] 
f.close() 
M = Matrix(data)

I then need to remove any particular scheduled slots that are not picked by any company:

M = matrix([row for row in M.transpose() if max(row) != 0]).transpose()

Once I've done this I can define the bi-partite graph (bi-partite simply means that the vertices can be separated in to 2 non adjacent collections):

g = BipartiteGraph(M)

We can then get a get a picture of this, I do this using a 'partition' (a graph colouring) that will colour the groups (red) and the schedule slots (blue):

g = BipartiteGraph(M)
p = g.coloring()
g.show(layout='circular',dist=1,vertex_size=250, graph_border=True, figsize=[15,15],partition=p)

The various options I pass to the `show` command are simply to get the circular arrangement (and other minor things):



The above looks quite messy and what I essentially want is get as many pairwise matchings between the blue vertices (slots) and red vertices (companies) so that each schedule slot is attributed at most 1 company and every company has at least 1 schedule slot.

On any given graph $G=(V,E)$ this problem is known as looking for a maximal matching and can be written down mathematically:

Max:  $\sum_{e \in E(G)} m_e$
Such that:  $\forall v$ $\sum_{e \in E(G) \atop v \sim e} m_e \leq 1$

We are in essence finding a subset of edges of our original graph in such a way as to maximise the number of edges such that no vertex has more than 1 edge.

This is all explained extremely well at the +Sage Mathematical Software System documentation pages here.

Furthermore at the documentation the code needed to solve the problem is also given:

p = MixedIntegerLinearProgram()  
matching = p.new_variable(binary=True)
p.set_objective(sum(matching[e] for e in g.edges(labels=False)))
for v in g:
      p.add_constraint(sum(matching[e]
          for e in g.edges_incident(v, labels=False)) <= 1)
p.solve()

When I run the above, `p` is now a solved Mixed Integer Linear Program (corresponding to the matching problem described). To obtain the solution:

matching = p.get_values(matching)
schedule = [e for e,b in matching.iteritems() if b == 1]

Calling `schedule` gives a set of edges (denoted by the corresponding vertex numbers):

[(5, 57), (0, 45), (23, 50), (4, 42), (38, 60), (26, 56), (34, 62), (16,
68), (1, 43), (7, 40), (9, 44), (36, 58), (12, 49), (35, 71), (28, 66),
(25, 47), (24, 53), (6, 46), (3, 64), (39, 67), (17, 69), (22, 55), (13,
48), (33, 41), (10, 63), (21, 61), (30, 52), (29, 65), (37, 70), (15,
54), (19, 51), (11, 59)]

It is then really easy to draw another graph:

B=Graph(schedule)
p = B.coloring()
B.show(layout='circular',dist=1,vertex_size=250, graph_border=True, figsize=[15,15],partition=p)

Which gives:



You can see that the obtained graph has all the required properties and most importantly is a lot less congested.

Some details.

I'm leaving out some details, for example I kept track of the names of the companies and also the slots so that the final output of all this looked like this:

4InARow: Fri1100
Abacus: Tue0930
Alpha1: Thu1130
Alpha2: Mon0930
AusfallLtd: Fri1230
AxiomEnterprise: Thu1000
Batduck: Thu1430
CharliesAngles: Thu1500
CwmniRhifau: Mon1330
EasyasPi: Fri1130
Evolve: Thu1300
HSPLLtd.: Mon1030
JECT: Tue1200
JJSL: Thu1030
JNTL: Mon1400
JennyCash: Tue1630
KADE: Fri1330
MIAS: Fri1300
MIPE: Thu1100
MLC: Tue1600
Nineties: Mon0900
Promis: Fri1400
R.A.C.H: Tue1530
RYLR: Tue1230
SBTP: Fri1030
Serendipity: Mon1230
UniMath: Tue1300
VectorEnterprises: Tue1330
Venus: Thu1400
codeX: Mon1300
dydx: Wed1630
eduMath: Thu0930

(BatDuck is my favourite company name by far...)

Why did I do this this way?

There are 3 reasons:

1. I shared the schedule with my students through a published Sage sheet on our server. That way they can see a direct applied piece of mathematics and can also understand some of the code if they wanted to.
2. "Point and click doesn't scale" - I'm sure I could have solved this one instance of my problem with pen and paper and some common sense faster than it took me to write the code to solve the problem. The thing is next year when I need to schedule these talks again it will at most take me 3 minutes as the code is all here and ready to go. (Most readers of this blog won't need that explained but if any of my students find their way here: that's an important message for you).
3. It was fun.