Linear Programming: can I formulate an objective to maximize multiple variables at once?
An answer to this question on Stack Overflow.
Question
Let's say I have some variables and constraints illustrated by the following system:
The gray lines can stretch and shrink an amount given by the range on top of them. The blue lines are just the endpoints and show how the gray lines interact.
My goal: I'd like to use linear programming to evenly maximize the size of the gray lines, like in the picture. You can imagine the gray lines with springs on them, all equally pushing outwards. A bad solution would be having all the blue lines pushed over as far to one side as possible. Note that there is a bit of leeway in this description, and multiple solutions are possible - all I need is for them to be reasonably even, and not have one value maxed out squishing everything else.
The objective function I tried simply maximizes the sum of line's size:
maximize: (B - A) + (C - B) + (C - A) + (D - C) + (E - B) + (E - D) + (F - E) + (F - D) + (F - A)
It's clear to me that this isn't a good solution, since the terms cancel out and an increase on one line just decreases it in another by the same amount, so the objective is never weighted towards evenly distributing the maximization among the variables.
I also tried to minimize each line's distance from their middle possible range. For line B - A, the middle value in it's range of (1,3) is 2. Here's the objective with the first term:
minimize: |(B - A) - 2| + ...
To implement the absolute value, I replaced the term with U and added additional constraints:
minimize: U + ...
with: U <= (B - A - 2)
U <= -(B - A - 2)
This has the same problem as the other objective: the difference is always proportional to the change in another line's difference. I think would work if I could square the difference, but I can't input that in linear solver.
Is there some objective function that would achieve what I am seeking, or is a linear solver just not the right tool for this?
I'm using Google OR-Tools, if that helps.
Here are the constraints written out:
1 <= B - A <= 3
0 <= C - B <= 1
1 <= C - A <= 4
9 <= E - B <= 11
7 <= D - C <= 11
0 <= E - D <= 1
3 <= F - E <= 5
3 <= F - D <= 6
15 <= F - A < = 15
Answer
Bear in mind that that your greatest problem is that you don't know what it is, exactly, that you want. So I've had to guess. Sometimes seeing a few guesses helps you refine what it is that you want, so this isn't too bad on your part, but it does make your question more difficult for the format of this site.
First, I'll assume that the springs can be modeled as a directed acyclic graph. That is, I can replace all the springs with arrows that point to the right. There will never be an arrow pointing from the right to the left (otherwise your springs would bend in a circle).
Once this is done, you can use set logic to figure out the identity of the leftmost blue bar. (I assume there is only one - it is left as an exercise to figure out how to generalize.) You can then anchor this bar at a suitable location. All the other bars will be positioned relative to it. This constraint looks like:
S[leftmost] = 0
Now, we need some constraints.
Each edge i has a source and end point (because the edges are directed). Call the position of the source point S and the position of the end point E. Further, the edge has a minimum length l and a maximum length L. Since we pin the location of the leftmost bluebar, the springs connected to it define the intervals in which their end points fall. Those end points are the source points for other springs, &c. Thus, each edge defines two constraints on the position of its end point.
S[i]+l[i] <= E[i]
E[i] <= S+L[i]
As an aside, note that we can now formulate a simple linear program:
min 1
s.t. S[leftmost] = 0
S[i]+l[i] <= E[i]
E[i] <= S+L[i]
If this program can be solved, then there is a feasible solution to your problem. Which is to say, the bar lengths don't produce a mutually inconsistent description of where the bluebars should be.
Now, we want to "evenly maximize the size of the gray lines", whatever this means.
Minimizing Deviation from Average Length
Here's one idea. The length the program chooses for each bar is given by E[i]-S[i]. Let's specify that this length should be "close to" the average length of the bar (L[i]+l[i])/2. Thus, the target quantity we want to minimize for each bar is:
(E[i]-S[i])-(L[i]+l[i])/2
Problematically, this value can be positive or negative depending on whether or not (E[i]-S[i])>(L[i]+l[i])/2. This isn't good because we want to minimize the deviation from (L[i]+l[i])/2, a value which should always be positive.
To cope, let's square the value and then take a square root, this gives:
sqrt(((E[i]-S[i])-(L[i]+l[i])/2)^2)
This might seem unsolveable, but stay with me.
Note that the foregoing is the same as taking the L2 norm of a one-element vector, so we can rewrite it as:
|(E[i]-S[i])-(L[i]+l[i])/2|_2
We can now sum the deviations for each bar:
|(E[0]-S[0])-(L[0]+l[0])/2|_2 + |(E[1]-S[1])-(L[1]+l[1])/2|_2 + ...
This gives us the following optimization problem:
min |(E[0]-S[0])-(L[0]+l[0])/2|_2 + |(E[1]-S[1])-(L[1]+l[1])/2|_2 + ...
s.t. S[leftmost] = 0
S[i]+l[i] <= E[i]
E[i] <= S+L[i]
This problem is not easily solveable in the form stated above, but we can perform a simple manipulation by introducing a variable t
min t[0] + t[1] + ...
s.t. S[leftmost] = 0
S[i]+l[i] <= E[i]
E[i] <= S+L[i]
|(E[i]-S[i])-(L[i]+l[i])/2|_2<=t[i]
This problem is exactly the same as the previous problem. So what have we gained?
Optimization is a game of converting problems into standard forms. Once we have a problem in a standard form, we can then Stand On The Shoulders Of Giants and use powerful tools to solve our problems.
The foregoing manipulation has turned the problem into a second-order cone problem (SOCP). Once in this form, it can be solved pretty much directly.
The code for doing so looks like this:
#!/usr/bin/env python3
import cvxpy as cp
import networkx as nx
import matplotlib.pyplot as plt
def FindTerminalPoints(springs):
starts = set([x[0] for x in springs.edges()])
ends = set([x[1] for x in springs.edges()])
return list(starts-ends), list(ends-starts)
springs = nx.DiGraph()
springs.add_edge('a', 'b', minlen= 1, maxlen= 3)
springs.add_edge('a', 'c', minlen= 1, maxlen= 4)
springs.add_edge('a', 'f', minlen=15, maxlen=15)
springs.add_edge('b', 'c', minlen= 0, maxlen= 1)
springs.add_edge('b', 'e', minlen= 9, maxlen=11)
springs.add_edge('c', 'd', minlen= 7, maxlen=11)
springs.add_edge('d', 'e', minlen= 0, maxlen= 1)
springs.add_edge('d', 'f', minlen= 3, maxlen= 6)
springs.add_edge('e', 'f', minlen= 3, maxlen= 5)
if not nx.is_directed_acyclic_graph(springs):
raise Exception("Springs must be a directed acyclic graph!")
starts, ends = FindTerminalPoints(springs)
if len(starts)!=1:
raise Exception("One unique start is needed!")
if len(ends)!=1:
raise Exception("One unique end is needed!")
start = starts[0]
end = ends[0]
#At this point we have what is essentially a directed acyclic graph beginning at
#`start` and ending at `end`
#Generate a variable for the position of each blue bar
bluevars = {n: cp.Variable(name=n) for n in springs.nodes()}
dvars = {e: cp.Variable() for e in springs.edges()}
#Anchor the leftmost blue bar to prevent pathological solutions
cons = [bluevars[start]==0]
for s,e in springs.edges():
print("Loading edge {0}-{1}".format(s,e))
sv = bluevars[s]
ev = bluevars[e]
edge = springs[s][e]
cons += [sv+edge['minlen']<=ev]
cons += [ev<=sv+edge['maxlen']]
cons += [cp.norm((ev-sv)-(edge['maxlen']-edge['minlen'])/2,2)<=dvars[(s,e)]]
obj = cp.Minimize(cp.sum(list(dvars.values())))
prob = cp.Problem(obj,cons)
val = prob.solve()
fig, ax = plt.subplots()
for var, val in bluevars.items():
print("{:10} = {:10}".format(var,val.value))
plt.plot([val.value,val.value],[0,3])
plt.show()
The results look like this:
If you want to "hand tune" the blue bars, you can modify the optimization problem we've built by adding weights w[i].
min w[0]*t[0] + w[1]*t[1] + ...
s.t. S[leftmost] = 0
S[i]+l[i] <= E[i]
E[i] <= S+L[i]
|(E[i]-S[i])-(L[i]+l[i])/2|_2<=t[i]
The larger w[i] is, the more important it will be that the spring in question is close to its average length.
Minimizing Squared Distance Between Ordered Blue Bars, Subject to Constraints
Using the same strategies as above, we can minimize the squared distance between the blue bars assume some sort of known order. This leads to:
min t[0] + t[1] + ...
s.t. S[leftmost] = 0
S[i]+l[i] <= E[i]
E[i] <= S+L[i]
|(S[i]-S[i+1])/2|_2<=t[i]
In the code below I first find feasible positions of the blue bars and then assume these map to a desirable order. Replacing this heuristic with more accurate information would be a good idea.
#!/usr/bin/env python3
import cvxpy as cp
import networkx as nx
import matplotlib.pyplot as plt
def FindTerminalPoints(springs):
starts = set([x[0] for x in springs.edges()])
ends = set([x[1] for x in springs.edges()])
return list(starts-ends), list(ends-starts)
springs = nx.DiGraph()
springs.add_edge('a', 'b', minlen= 1, maxlen= 3)
springs.add_edge('a', 'c', minlen= 1, maxlen= 4)
springs.add_edge('a', 'f', minlen=15, maxlen=15)
springs.add_edge('b', 'c', minlen= 0, maxlen= 1)
springs.add_edge('b', 'e', minlen= 9, maxlen=11)
springs.add_edge('c', 'd', minlen= 7, maxlen=11)
springs.add_edge('d', 'e', minlen= 0, maxlen= 1)
springs.add_edge('d', 'f', minlen= 3, maxlen= 6)
springs.add_edge('e', 'f', minlen= 3, maxlen= 5)
if not nx.is_directed_acyclic_graph(springs):
raise Exception("Springs must be a directed acyclic graph!")
starts, ends = FindTerminalPoints(springs)
if len(starts)!=1:
raise Exception("One unique start is needed!")
if len(ends)!=1:
raise Exception("One unique end is needed!")
start = starts[0]
end = ends[0]
#At this point we have what is essentially a directed acyclic graph beginning at
#`start` and ending at `end`
#Generate a variable for the position of each blue bar
bluevars = {n: cp.Variable(name=n) for n in springs.nodes()}
#Anchor the leftmost blue bar to prevent pathological solutions
cons = [bluevars[start]==0]
#Constraint each blue bar to its range
for s,e in springs.edges():
print("Loading edge {0}-{1}".format(s,e))
sv = bluevars[s]
ev = bluevars[e]
edge = springs[s][e]
cons += [sv+edge['minlen']<=ev]
cons += [ev<=sv+edge['maxlen']]
#Find feasible locations for the blue bars. This is a heuristic for getting a
#sorted order for the bars
obj = cp.Minimize(1)
prob = cp.Problem(obj,cons)
prob.solve()
#Now that we have a sorted order, we modify the objective to minimize the
#squared distance between the ordered bars
bar_locs = list(bluevars.values())
bar_locs.sort(key=lambda x: x.value)
dvars = [cp.Variable() for n in range(len(springs.nodes())-1)]
for i in range(len(bar_locs)-1):
cons += [cp.norm(bar_locs[i]-bar_locs[i+1],2)<=dvars[i]]
obj = cp.Minimize(cp.sum(dvars))
prob = cp.Problem(obj,cons)
val = prob.solve()
fig, ax = plt.subplots()
for var, val in bluevars.items():
print("{:10} = {:10}".format(var,val.value))
plt.plot([val.value,val.value],[0,3])
plt.show()
That looks like this:
Minimizing Squared Distance Between All Blue Bars, Subject to Constraints
We could also try to minimize all of the pairwise squared distances between blue bars. To my eye this seems to give the best result.
min t[i,j] + ... for all i,j
s.t. S[leftmost] = 0
S[i]+l[i] <= E[i] for all i
E[i] <= S+L[i] for all i
|(S[i]-S[j])/2|_2 <= t[i,j] for all i,j
That would look like this:
#!/usr/bin/env python3
import cvxpy as cp
import networkx as nx
import matplotlib.pyplot as plt
import itertools
def FindTerminalPoints(springs):
starts = set([x[0] for x in springs.edges()])
ends = set([x[1] for x in springs.edges()])
return list(starts-ends), list(ends-starts)
springs = nx.DiGraph()
springs.add_edge('a', 'b', minlen= 1, maxlen= 3)
springs.add_edge('a', 'c', minlen= 1, maxlen= 4)
springs.add_edge('a', 'f', minlen=15, maxlen=15)
springs.add_edge('b', 'c', minlen= 0, maxlen= 1)
springs.add_edge('b', 'e', minlen= 9, maxlen=11)
springs.add_edge('c', 'd', minlen= 7, maxlen=11)
springs.add_edge('d', 'e', minlen= 0, maxlen= 1)
springs.add_edge('d', 'f', minlen= 3, maxlen= 6)
springs.add_edge('e', 'f', minlen= 3, maxlen= 5)
if not nx.is_directed_acyclic_graph(springs):
raise Exception("Springs must be a directed acyclic graph!")
starts, ends = FindTerminalPoints(springs)
if len(starts)!=1:
raise Exception("One unique start is needed!")
if len(ends)!=1:
raise Exception("One unique end is needed!")
start = starts[0]
end = ends[0]
#At this point we have what is essentially a directed acyclic graph beginning at
#`start` and ending at `end`
#Generate a variable for the position of each blue bar
bluevars = {n: cp.Variable(name=n) for n in springs.nodes()}
#Anchor the leftmost blue bar to prevent pathological solutions
cons = [bluevars[start]==0]
#Constraint each blue bar to its range
for s,e in springs.edges():
print("Loading edge {0}-{1}".format(s,e))
sv = bluevars[s]
ev = bluevars[e]
edge = springs[s][e]
cons += [sv+edge['minlen']<=ev]
cons += [ev<=sv+edge['maxlen']]
dist_combos = list(itertools.combinations(springs.nodes(), 2))
dvars = {(na,nb):cp.Variable() for na,nb in dist_combos}
distcons = []
for na,nb in dist_combos:
distcons += [cp.norm(bluevars[na]-bluevars[nb],2)<=dvars[(na,nb)]]
cons += distcons
#Find feasible locations for the blue bars. This is a heuristic for getting a
#sorted order for the bars
obj = cp.Minimize(cp.sum(list(dvars.values())))
prob = cp.Problem(obj,cons)
val = prob.solve()
fig, ax = plt.subplots()
for var, val in bluevars.items():
print("{:10} = {:10}".format(var,val.value))
plt.plot([val.value,val.value],[0,3])
plt.show()
That looks like this:
