Spring chain example

Let’s say that you want to compute the balance positions of a chain of springs.

Start with includes:

#include <iostream>
#include <QuadProgMm.hpp>

using namespace QuadProgMm;

Define a class representing a spring. It’s a simple aggregate with the spring’s strength, and unconstrained length:

struct Spring {
  Spring(float strength_, float length_) :
    strength(strength_), length(length_)
  {}

  float strength, length;
};

Then, a class for the chain. It has two boundaries (left and right), a vector of Spring, and a vector of Variable for the positions to be computed:

struct SpringChain {
  SpringChain(
    float left_, float right_,
    const std::vector<Spring>& springs_
  ) :
    left(left_), right(right_),
    springs(springs_),
    positions(makePositions(springs_)),
    solution(resolve())
  {}

  double get(size_t i) const {
    return solution.get(positions[i]);
  }

  private:
    static std::vector<Variable> makePositions(
      const std::vector<Spring>&
    );
    Solution resolve();

    float left, right;
    std::vector<Spring> springs;
    std::vector<Variable> positions;
    Solution solution;
};

Then, the definition of makePositions that constructs an independent Variable for each position to solve:

std::vector<Variable> SpringChain::makePositions(
  const std::vector<Spring>& springs
) {
  std::vector<Variable> positions;
  for(size_t i = 0; i != springs.size() + 1; ++i) {
    positions.push_back(Variable());
  }
  return positions;
}

Then, the most important part: the definition of SpringChain::resolve:

Solution SpringChain::resolve() {

First, define the objective: each spring has a potential energy of \(1/2 \cdot k \cdot (l - l_0) ^ 2\), where \(k\) is its strength, \(l_0\) is its unconstrained length, and \(l\) is its current length. The balance position is reached when the total energy is minimal: the objective is to minimize the sum of all those energies. So, build the total energy as a quadratic expression of the variables:

QuadraticExpression energy = 0;
for(size_t i = 0; i != springs.size(); ++i) {
  LinearExpression l_minus_l0 =
    positions[i + 1] - positions[i]
    - springs[i].length;

  energy += 0.5 * springs[i].strength * l_minus_l0 * l_minus_l0;
}

Add the constraints on boundaries:

std::vector<Constraint> boundaries {
  positions.front() == left,
  positions.back() == right,
};

And solve the QP problem:

  return minimize(energy, boundaries);
}

Finally, use the SpringChain class:

int main() {
  std::vector<Spring> springs {
    Spring(1, 2),
    Spring(1, 3),
    Spring(10, 5),
    Spring(1, 2),
  };

  SpringChain chain(0, 10, springs);
  for(size_t i = 0; i != springs.size() + 1; ++i) {
    std::cout << "Position " << i << ": "
      << chain.get(i) << std::endl;
  }
}

And here is the result:

Position 0: 7.60026e-10
Position 1: 1.35484
Position 2: 3.70968
Position 3: 8.64516
Position 4: 10