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