Computing forces in a system of beams

*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:

- The sum of the forces acting on each beam is zero.
- The sum of the torques of the forces acting on each beam is zero.
- The sum of the forces acting on each hinge is zero.

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

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:

- The hinge connects the beam to the ground. In this case, there is an extra unknown ground force, which we don't care about, so this connection doesn't produce a useful equation. We could add both the unknown force and the balance equation to the system, which would only increase the computational efforts without adding any useful information.
- A beam end doesn't connect to anything and simply floats in space. In this case it isn't a hinge or a connection of any sort, but it is still useful to consider it as a single-beam hinge, and the hinge equation would tell us that the force at this beam's free end must be zero (because it is static and isn't connected to anything).

So, here's our final mathematical model of a set of beams:

- For each beam, we have 2 unknown 2D vectors \( F_a, F_b \) representing the forces acting on the beam's ends
- For each beam, we have a vector equation \( F_a + F_b + mg = 0 \)
- For each beam, we have a scalar equation \( F_a \cdot n = F_b \cdot n \)
- For each hinge, we have a vector equation \( \sum F_i = 0 \) (
*for this equation, wall connections aren't hinges, while free-floating beam ends are hinges*)

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:

- For each beam, we have 2 unknown 2D vectors, thus 4 scalar unknowns
- For each beam, we have 1 vector equation and 1 scalar equation, so 3 scalar equations in total
- For each hinge, we have 1 vector equation, which is 2 scalar equations

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:

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:

- Improve our model to better reflect reality, so that it always produces a non-singular square system of equations
- Solve our system in a least-squares sense, finding an arrangement of forces that is as close as possible to satisfying the equations
- Slap a stupid iterative equation solver and hope for the best

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:

- Take an initial guess of the solution (setting all unknowns to zero works fine)
- Iterate over all equations you have
- For each equation of the form \( X+Y+Z=A \), compute the
*current*value \( A' = X+Y+Z \), compute the*error*\( A - A' \), and change the values of the unknowns by evenly redistributing the error \( X \leftarrow X + (A-A')/3 \) (and so on for \( Y \) and \( Z \)), forcing this particular equation to be true - Repeate 2-3 until your unknowns converge to a solution

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:

- For overdetermined systems with no proper solution, it produces some set of forces, with no guarantee on how good this forces are
- For underdetermined systems it produces some solution which might look pretty random if the initial guess is bad, or if we're adding beams on-the-fly

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 :)