1

Let's say I have some variables and constraints illustrated by the following system: enter image description here 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
Richard
  • 56,349
  • 34
  • 180
  • 251
lifeformed
  • 525
  • 4
  • 16
  • It seems you did not take your time to process the answers to your earlier question. Then, probably you will recognize that your objective-description is either very informal or just conflicting (`evenly maximize the size of the gray lines` vs. `The objective function I am using tries to maximize each line's size`) – sascha Oct 25 '18 at 11:10
  • You say, you're using Google or-tools. Where is the code then? – mattmilten Oct 25 '18 at 16:36
  • @sascha That other question was about a fast way of finding the bounds of the individual terms; this one is a new question about finding a good solution to the system as a whole. For this, my goal is to evenly maximize the sizes of each line. That objective function is just what I thought would achieve that goal. I will clarify in an edit. Thanks for your response to my other question by the way, I did indeed process it and it helped me come to my solution. – lifeformed Oct 25 '18 at 18:48
  • @mattmilten I edited the post with the details formalized a little. I didn't include the code because it's just boilerplate on top of my description, and my question is more about the concept than specifics of the code. I only mentioned OR-tools in case maybe it has some capabilities I overlooked. – lifeformed Oct 25 '18 at 18:49
  • `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` I don't get this. Reducing A-B by 1 has the potential to gain 1 + 1 with B-E and B-C, simply because of the asymmetric cardinality. – sascha Oct 25 '18 at 18:52
  • That's true, I didn't consider that. However, it would still give an undesirable result, since it prioritizes connectedness over evenness. In this example, `D-C` would be squished to it's minimum since it's less connected than it's neighbors (which actually wouldn't be too bad in this example, but it could in others). It's okay if a term is punished for low connectivity, but I don't want it to be completely minimized, rather only shortened proportional to the connectivity. Also, when the number of neighbors on either side is even, then it wouldn't consider even spacing at all. – lifeformed Oct 25 '18 at 19:26
  • I don't have a good sense of what you're trying to optimize or what the number in your picture mean. – Richard Oct 26 '18 at 20:11
  • @Richard Imagine the gray lines are metal bars that can shrink and expand. The numbers on them represent the minimum and maximum size they can shrink and expand to. The vertical blue lines are plates that connect the bars together. The solution I want to find is a state where the bars have expanded at the same time and settled in equilibrium. Imagine there are springs around all the bars, pushing outwards. I want the final resting state of that spring-based system. – lifeformed Oct 26 '18 at 20:46

1 Answers1

5

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:

Location of terminii in spring optimization problem

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:

Ordered lines in optimization problem

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:

Unordered lines in optimization problem

Community
  • 1
  • 1
Richard
  • 56,349
  • 34
  • 180
  • 251
  • 1
    Thanks so much for your detailed reply. I am going to experiment with the code you posted and get back to you tomorrow. You are correct in your understanding of my goal, by the way. My first thoughts are this: I attempted a similar solution, but I ran into a problem: minimizing the difference from the midpoint would sometimes give results where one bar was "sacrificed", squished to its minimum width, while the other bars would then have more leeway to get close to the midpoint. (cont.) – lifeformed Oct 27 '18 at 08:52
  • That solution was just as valid as one where they all equally tried to get to the midpoint but were each off by a smaller amount. It's most readily apparent with this simple system: A-----B-----C-----D, with each range being between 1 to 13, and the total distance (D-A) being fixed to 15. The ideal solution is A-----B-----C-----D, (distance of 5 between all terms), but A-------B-------C-D (distance of 7, 7, and 1 respectively) is equally valid: the sum of their differences from 7 (the midpoint of 1 to 13) is 6 in both solutions. The solver I used would return the bad solution more often. – lifeformed Oct 27 '18 at 09:05
  • @lifeformed: It is not general, but in these cases you could try adjusting the weights or adding a regularisation term. – Richard Oct 27 '18 at 16:17
  • @lifeformed: I've taken a couple more stabs at it. But if you're getting bad results on some of your test cases, it's probably good if you upload the values for one of those. However, the formulations I've provided, when weighted (see answer), provide a fairly general framework. – Richard Oct 28 '18 at 02:02
  • Thanks! One problem I am encountering is that the linear solver in Google OR-Tools doesn't appear have a way to do squares as a constraint. Is a good trick to approximate it? I tried adding more linear constraints that kind of give a square-like curve but I'm not sure if that's good practice. – lifeformed Oct 29 '18 at 03:24
  • @lifeformed: I couldn't find documentation demonstrating that Google OR supports SOCPs. I think it might be restricted to linear programs, in which case you'd have to find a way to approximate the square. I'd suggest using `cvxpy` instead: it uses cvxopt under the hood, which can do both LPs and SOCPs. – Richard Oct 30 '18 at 15:48
  • 1
    Thanks for your help, this has set me down the right path. – lifeformed Nov 03 '18 at 07:39