What I'm describing in this post turned out to be a bad mathematical model, since I've missed one crucial component. If you want to skip this and learn about the actually good model, go straight to this article.
There is one idea I've been passively tinkering with in my head for years without any success. Imagine a game where you have some building blocks (cubes, walls, platforms, pillars, whatever), and you want them to be able to break under load. Imagine Minecraft, where you build a horizontal beam of blocks in the air, attached to some mountain, and the beam breaks if you build it too far. Or, say, the Sims, where a small wooden pillar breaks because your whole 5-floor concrete house was supported by this only pillar. You get the idea.
Game design issues aside, how do you even calculate this load on some supporting structure? We could probably do something simple, like if you have 5 pillars connecting the house to the ground, and the house weighs 50 tons, the load is 10 tons per pillar. This, however, ignores a lot of detail: what about supports that are not connected to the ground and instead connect different floors of the house? What about uneven mass distribution across the house? What about the center of gravity and how it affects the distribution of load? Etc, etc.
After thinking about Minecraft-like blocks for a while, I realized this is already too complicated, and I don't understand a lot of things involved here. (Though, I've managed to come up with a really nice model, which I'll probably describe in a follow-up post.) Then I decided to constrain myself to a hopefully much simpler case of beams.
By the way, of course I tried to research the topic from some physics books and didn't understand anything. I also tried to ask something related to this on physics.stackexchange, got three answers, and didn't understand a single one of them either.
So, I have some configuration of beams. They can connect to walls, or the floor, or to each other. Something like this:
Let's say they are all uniformly dense, and have the same density. We treat them as one-dimensional, so the density is in units mass per length. They are all affected by gravity.
Whenever a beam touches a wall/floor, or whenever two beams connect to each other via their ends, let's say this connection is done via an incredibly strong glue, so that the beams cannot rotate or somehow move around this connection. This connections are rigid and massless. I'm calling them hinges in my code, and despite them not being hinges, I'll still call them hinges throughout this post.
We also assume that the beams are static – they don't move, and are in an equillibrium (even if it is an unstable one).
What we want to compute are the forces acting on each beam's two ends. Something like this:
From that, as we'll see later, we can compute how the beam is stretched or contracted, how much is it bent, and so on.
You might say: there is nothing simpler! Just write down the equations for these forces, solve them, and you're done! Well, you'd be mostly right, except for the solving part.
But before we discuss why is it so hard to solve these equations, let's look at some examples first.
Say we have a single perfectly straight vertical beam. Its lower end touches the ground, while its upper end doesn't touch anything. The force on the upper end must be equal to zero, because there is nothing at this end acting on the beam (we ignore things like air pressure, wind, etc). Thus, the force at the lower end must point upwards and completely compensate the gravity force:
A similar example, but this time the beam is glued to the ceiling. The bottom force must be zero, and the top force must again point upwards to compensate for gravity:
Another, somewhat silly, example: a beam is simply lying on the floor. However, in our model it can only glue to the floor on its ends, while the middle section of the beams is kind of above the floor by an infinitely small amount. In this case, the gravity must be compensated by two forces on the left and right ends, which should probably be equal due to symmetry:
Ok fine, but what if our beam is sitting diagonally between a floor and a wall? Honestly, I don't have a good intuition about what should happen in this case. The forces must compensate the gravity, but other than that I'm pretty lost concerning how they should look like. Here's a fairly reasonable arrangement of forces:
Note that the forces don't have to be orthogonal to the walls, as would happen if the beam was touching the wall without any glue.
What about more than one beam? Say we have two beams arranged in an upside-down V shape. They must be pushing each other sideways at the connection point, and they probably push the ground sideways as well, and the ground pushes them back. Something like this:
Or consider this nice symmetric triangular arrangement of beams, kind of like a triangular house of cards:
Or maybe we go crazy and make a whole bunch of beams all touching the ground and connected to a single point:
I hope you have an idea of what we're trying to calculate by now. Let's try making our model precise now.
So, we want to compute the forces at both ends for each beam. This means that for each beam we have two unknown 2-dimensional force vectors \( F_a \) and \( F_b \).
We need some equations for these forces, which come from the assumption that our system is in a static equillibrium. We have three types of equations:
There are only three forces acting on each beam: the two forces \( F_a, F_b \) acting on the beam's ends, and the gravity force \( mg \). So, the first equation tells us that
\[ F_a + F_b + mg = 0 \]The torque, i.e. rotational moment, of a force \( F \) is computed as \( r \times F \), where \( r \) is the vector from the center of the body to the point where the force is applied to, and \( \times \) is a cross product. The gravitational force never does any rotation (this is actually a nontrivial result from mechanics; we can think that gravity always acts on the body's center, in which case \( r = 0 \)), so we're only left with the forces acting on the beam ends. Since the ends are symmetric with respect to the center of the beam, the \( r \) vectors are opposite for \( F_a \) and \( F_b \), and we get something like
\[ F_a \times (-r) + F_b \times r = 0 \]or
\[ (F_b - F_a) \times r = 0 \]Now, in 2D the cross product can be replaced with a dot product with a vector, orthogonal to \( r \). In our case \( r \) is the vector from the beam center to the beam's end, so an orthogonal vector is, for example, a vector \( n \) which is normal to the beam. We don't care about it's length, because our equation has a \( = 0\) in the right-hand side. So, our equation becomes
\[ (F_b - F_a) \cdot n = 0 \]or
\[ F_a \cdot n = F_b \cdot n \]which means that in order not to rotate the beam, the two forces must coincide in the normal direction. This makes total sense: if one force was stronger than the other, it would overcome it and rotate the beam. This also has the added benefit that the torque balance equation has the dimension of forces instead of forces times lengths, which will be useful later.
The last equation is the sum of all the forces for each hinge. When a set of beams are connected in a single point, the sum of the forces acting this point must be zero. For example, if three beams are connected to a single hinge, with the first two being connected via their first end (where \( F_a \) acts), and the third being connected via its second end (where \( F_b \) acts), this gives us the equation
\[ F_{a1} + F_{a2} + F_{b3} = 0 \]where \( F_{a1} \) is the \( F_a \) force of the first beam, and so on.
I've glossed over two special cases:
So, here's our final mathematical model of a set of beams:
Now let's talk about how we'd solve this equation.
The first thing that comes to mind about this set of equations is that it looks stupidly simple. A few sums here and there are equal to zero, that's it. Don't let these equations fool you: they are cursed as hell. Well, maybe not that cursed, but still far from being completely tame.
The second thing that comes to mind about this set of equations is that it is linear. The components of our unknown forces are only added up and multiplied by scalars, that's it. This is insanely good: linear systems are basically the only thing humanity knowns how to reliably solve. If the system is good enough, that is. Is it?
For a linear system to be good, it must be square: the number of unknowns must be equal to the number of equations. Let's count how many of these we have:
In total, if we have \( B \) beams and \( H \) hinges, we have \( 4B \) unknowns and \( 3B+2H \) equations. Solving \( 4B=3B+2H \) gives us \( B=2H \), i.e. there must be twice as many beams as there are hinges. Must this always be the case? Hardly so!
Look at our first example, a beam standing vertically on the ground.
We have \( B=1 \) and \( H=1 \) (the free top end of the beam is a virtual hinge, remember?). In this case, we get \( 4B = 4 \) unknowns and \( 3B+2H=5 \) equations. Tsk-tsk.
Consider another example of a diagonal beam supported at both ends.
\( B=1 \) beam and no hinges (the wall hinges don't lead to useful equations), so we get \( 4B=4 \) unknowns and only \( 3B+2H=3 \) equations.
Ok, but can we get a square system? Yes, sometimes, like in the upside-down V shape:
Here, \( B=2 \) and \( H=1 \), leading to \( 4B=8 \) unknowns and \( 3B+2H=6+2=8 \) equations. Hooray!
Why is the number of unknowns and equations so important, you might ask? Well, because it tells us how determined our system is.
I'm not even sure if "determinacy" is a word, but I like it.
Update: I googled it. It is, in fact, a word.
If the system is square, i.e. the number of equations is the same as the number of unknowns, basic linear algebra tells us that there is usually a single solution to the system. (We'll talk about this "usually" in a bit.)
If the number of equations is smaller than the number of unknowns, we have an underdetermined system: we didn't constrain our unknowns enough, and they have some wiggle room. In this case, usually there are many different solutions. Infinitely many, even.
If, however, the number of equations is larger than the number of unknowns, the system is overdetermined: we have more constraints than we can possibly satisfy. In this case, the system usually doesn't have a solution at all. Like, at all. No solution, period. Well, usually.
This is all somewhat shady, because we're talking about a physical system. Surely, it is a model of a system, but a seemingly good model nevertheless (if you know a better model of the same thing – email me immediately!). Surely the physical world can compute some forces, but for some reason we can't? That is pants, if you ask me.
The case of an underdetermined system is easy to justify: probably the exact force distribution depends on the specific materials, on some friction coefficients, and so on. We don't include them in our model, hence there is no unique solution to the problem.
The case of an overdetermined system is harder. I mean, I have no idea why a system like this would be overdetermined, yet most interesting systems seem to be of that type. Somehow we included more equations that physics itself allows for, and I'm struggling to find a reasonable interpretation for that, except that maybe our model is just bad and wrong.
Or maybe I'm just not that good at physics.
Spoiler: I am indeed bad at physics, and as I've said in the beginning, the model is indeed bad and wrong, and I've already figured a better one.
What's up with that usually in the previous section? Well, even if a system is square, it might still be bad if the determinant of the matrix of the system is zero, or, equivalently, if the equations are not independent. In this case, we call the system singular, and the system can have no solutions (like for an overdetermined system), or have many solutions (like in an underdetermined system). Here's an example of how this might happen:
\[ \begin{matrix} x+y=1 \\ 2x+2y=2 \end{matrix} \]Clearly, the second equation follows from the first one, so we could as well just throw it away and be left with 1 equation for two unknowns, which has infinitely many solutions.
For another example, consider a system
\[ \begin{matrix} x+y=1 \\ 2x+2y=3 \end{matrix} \]Clearly, the second equation contradicts the first one, and there can be no solutions.
This only happens when the system has zero determinant, though, and a random system in the wild typically has a non-zero determinant.
However, this can also happen to an underdetermined or an overdetermined system, in which case the determinant doesn't make sense and we talk about the system having or not having a full rank.
Our systems come from a set of beams, and whether they are full-rank or not depends on a particular system, i.e. a particular arrangement of beams.
So, we have a system of linear equations, which might or might not have a solution, and which might have one or infinitely many solutions. Great.
What we want is to solve it, i.e. to find some arrangement of forces that satisfies all the equations, or at least satisfies them as much as possible. We have a few options to achieve this:
Guess what I tried first.
There's a really funny way of solving linear systems which works far more often than you'd think it should. Here's the recipe:
As far as I've heard, this is called iterative relaxation or something like that. In some cases, like solving a Laplace's equation on a square grid, it directly corresponds to iterations of the Gauss-Seidel method (or at least I think it does...), so at least this method isn't completely hopeless.
This solver doesn't take any assumptions about our system of equations, and is extremely easy to implement. The only problem is the torque balance equation, which is of the form
\[ (F_b-F_a)\cdot n = 0 \]instead of being a simple sum. In this case, I'm computing the current value \( M = (F_b - F_a)\cdot n \), and alter the forces like this:
\[ \begin{matrix} F_a \leftarrow F_a + \frac{M}{2} n \\ F_b \leftarrow F_b - \frac{M}{2} n \end{matrix} \]If we compute the new value, we get
\[ \left[\left(F_b-\frac{M}{2}n\right) - \left(F_a+\frac{M}{2}n\right)\right]\cdot n = (F_b-F_a)\cdot n - 2\frac{M}{2}(n\cdot n) = M - M = 0 \](We used the fact that \( n \) is normalized).
So, this is an iterative method, and on each step it tries to force each equation to be true, one by one. Let's see how well it works.
First, a square system: it seems to converge to a reasonable solution pretty quickly:
Let's look at a few overdetermined systems – one that has a solution (a straight vertical beam), and one that doesn't (a beam that is diagonal but is still supported only by the floor):
The method manages to find a solution when it exists, even though the system is overdetermined, and also finds some solution if no proper solution is possible.
Let's have a look at an underdetermined system: a beam supported both by the floor and the wall. In fact, let's have two such beams simultaneously.
The method seems to find a good enough solution! But, by the nature of the method, if we start from a different initial force arrangement, it will converge to a different solution. Here, I've manually set the starting forces to something non-zero:
This example is interactive, by the way! Click and drag to add new beams, right-click to remove beams. You can try to recreate the exact same beams and you'll see that the solution will depend on how exactly did you create this beams.
Let's also stress-test this method:
It does something reasonable (again, provided we zero-initialize all the forces), but it takes some time to converge.
I was still generally happy with this method, except for how it handles non-square systems:
Can we do something better?
There's a well-known trick of turning any equation into a minimization problem: instead of solving \( f(x)=0 \), you minimize \( f(x)^2 \). We square the thing so that it becomes positive and \( f(x)=0 \) will be the true minimum, and we don't do something like \( \|f(x)\| \) because squaring is a simple algebraic operation with nice properties. The benefit is that even if the original problem didn't have a solution, the minimization problem almost always has some solution.
So, in our case, instead of solving the overdetermined linear equation \( Ax=b \) which might not even have a solution, we try to minimize \( \|Ax-b\|^2 \). Note that this is a sum of squares of all the coordinates of the vector \( Ax-b \), and if these coordinates have different units, adding them up doesn't quite make sense. Here, though, we've already transformed the torque balance equations to have the same units as the force balance equations.
There are explicit formulas for that, which, however, require us to solve a square linear system in the process (using e.g. Gauss elimination), and I wanted to implement something in like 15 minutes to see if the whole idea can work at all. So, instead I decided to try gradient descent.
If you stare long enough into the formula \( \|Ax-b\|^2 \), you can figure out that the gradient of this scalar function with respect to vector \( x \) is equal to \( 2 A^T (Ax-b) \). Thus, we can try doing the usual gradient descent with a fixed step size (I used 0.25) and see what solution we get. This seems to work better (right image) for overdetermined systems than the stupid solver (left image):
This method found a better solution (i.e. the total error \( \|Ax-b\|^2 \) was smaller) at the expense of ignoring the physical reality and assigning some force to the free-floating top beam end.
For square systems, this method performs as good as the previous one, though it seems to converge a bit slower:
And for underdetermined systems this method is still sensitive to the initial guess or dynamic modifications of the system:
(This simulation is also interactive.)
And, as usual, a stress-test:
It converges, but slower than the previous method.
We still haven't solved the problem with underdetermined systems. It would be cool if we could just minimize something by another gradient descent, and I've seen in John D. Cook's blog that you can formulate an underdetermined least-squares problem as minimizing \( \|Ax-b\|^2 + \|x\|^2 \). It kind of makes sense at first glance: for a solution to \( Ax=b \) the \( \|Ax-b\|^2 \) term will be zero, and we'll find a minimum-length solution by minimizing the \( \|x\|^2 \) term. The gradient changes only slightly compared to the previous section: \( 2A^T (Ax-b)+2x \), and if we try to apply it to an underdetermined system, we get
As you can see, the error \( \|Ax-b\|^2 \) stays positive, i.e. the method converges to something that is not a solution to \( Ax=b \), even though the original problem is underdetermined and has infinitely many solutions!
When I asked about it on twitter, people got confused, thinking I'm asking about some ML-related regularization. A bit later we figured out that minimizing \( \|Ax-b\|^2+\|x\|^2 \) is not equivalent to solving \( Ax=b \), but to some different equation:
That's what you get for relying on random blogs that don't cite sources, I guess.
Yes, I don't cite sources either. Oh well.
So, I abandoned the idea of having a simple iterative solver for the underdetermined case, and instead decided to do thing the right way, which means solving the normal equations for our least-squares problem.
For the overdetermined case, minimizing \( \|Ax-b\|^2 \) explicitly leads to the formula
\[ x = (A^T A)^{-1} A^T b \]For the underdetermined case, we need to reformulate the problem as minimizing \( \|x\|^2 \) under the constraint \( Ax=b \), which sounds different but leads to a very similar (but still different!) formula
\[ x = A^T (A A^T)^{-1} b \]Both these methods require inverting a matrix or, alternatively, solving a square linear system. If our system \( Ax=b \) is square from the start, we just apply Gauss elimination to it.
This is not an iterative method, so we don't get the forces smoothly converging to a solution, but instead we get the solution instantly. Also, nor does Gauss elimination exploit the fact that our matrix \( A \) is very sparse, – most forces only participate in a few equations, – and neither does it care about \( AA^T \) or \( A^TA \) being symmetric positive-semidefinite. To exploit these, we could do something like sparse Cholesky decomposition of \( AA^T \) or \( A^TA \) and then use it to perform elimination to solve our system.
So, a few sparse matrix multiplications later, I got this method to work, and it seems to be pretty good. For the overdetermined system with no solution, it computes a solution as good as the gradient descent:
and for an underdetermined system, it tends to compute the solution which is the most symmetric and reasonable:
In fact, all of the systems in the examples section were computed using this method, so I won't show too many examples here. Have a look at this majestic bridge, though:
By the way, if we remove a few beams from this bridge, we get a singular system – something we couldn't detect with our previous method, but the Gauss elimination can easily report that:
If you remember the beginning of this post, I actually wanted to calculate the stresses on each beam, to figure out if it breaks or not. I will define bending stress and contraction/stretching stress.
Bending is how much do the forces acting on the beam, well, bend it. The stronger the forces in the normal direction, the higher the bending. Since we know that \( F_a\cdot n = F_b \cdot n \), we can define either of these values as the bending stress \( B \).
Now let's add \( F_a\cdot n \) and \( F_b\cdot n \). They are equal, thus the sum is equal to \( 2B \). However, we also have the force balance equation \( F_a+F_b+mg=0 \), leading to \[ 2B = (F_a+F_b)\cdot n = -mg \cdot n \]
Thus, the amount of bending is completely determined by the gravity force and the beam's geometry (the vector \( n \)), and doesn't depend on the forces we've been trying so hard to compute. Therefore, let's pretend we're not interested in bending that much :)
Now let's look at stretching/contraction stresses. These only depend on the tangential components of the forces, i.e. the components that go parallel to the beam (as opposed to the bending depending on the components that are orthogonal to the beam). If the end force \( F_b \) points towards the beam center and is opposite to the beam direction vector \( t \), while the start force \( F_a \) points again towards the beam center and is aligned with the beam direction vector \( t \), we get contraction: two forces acting on the beam from both sides are trying to squeeze it.
If, however, both \( F_a \) and \( F_b \) point outwards, we have stretching: the forces are trying to stretch the beam.
We can figure out which one is the case simply by computing \( (F_b - F_a) \cdot t \): if this number is positive, we have stretching, otherwise we get contraction.
Here I'm showing stretching/contraction using inward/outward white arrows:
Here's an interactive playground where you can try different beam configurations. Left click and drag to create beams (they automatically connect to each other via hinges if their endpoints coincide), right click to remove beams. You can select whether to show forces or stresses using the control panel below the thing, and also change the scale of the force vectors. Don't hesitate to show me some funny configurations (in twitter, mastodon, or simply be email), I'd love to see what you come up with!
| Force scale: |
So, here's our model for computing the forces in a beam. I'm not quite satisfied with it, but it seems to work in most cases. Maybe a more physically appropriate way of doing this would involve minimizing something smarter than just the sum of squared errors.
I guess what most bridge-building and similar games do is actually simulate the physics (they have dynamic physics, after all!), and look at the constraint forces that the physics engine computed. It would be interesting to see how these would differ from my naive least-squares method :)
Anyway, I hope you liked the post. I'm planning to write another one, experimenting with squares or cubes this time. Stay tuned :)