Archive | Software RSS for this section

Finding a Better Truss

This was an experiment in using simulated annealing to improve the load bearing capacity of a truss.

Essentially, SA is a technique used to find the “best” state of a given system. It stems from the metallurgical technique of heating and cooling a metal to reduce the Gibbs energy and relieve internal stresses. In general terms we begin with a problem which has a certain energy landscape, and use SA to find either the highest or lowest point on the landscape. The animation below (Wikipedia) is probably the simplest explanation:

Simulated Annealing

In Layman’s Terms

If you start at a point on the landscape, say somewhere in the leftmost valley above, you can easily find a “better” spot by just walking up the hill. But eventually you’ll get to the top and be unable to move any higher, even though you can clearly see that gigantic mountain in the middle. To get there, you have to go back down! But at least you can see the mountain peak so you know you’re not in the best spot, because the best spot is over there. What if you were blind, and all you had was an altimeter to tell you how high you were? You wouldn’t be able to see the mountain, so from your perspective you would have already found the best spot. SA is just what you need to find the actual best spot.

Our blind mountain climber with an altimeter – let’s call him Dave – is a smart guy. He knows that if his altimeter shows a certain value, and whenever he moves it shows a lower value, he’s at a high point: a local maximum. But he wants to find the single highest point: the global maximum. Dave knows that he will sometimes have to move down before he can start moving back up, and that when he does he may end up in a better position.

At dawn, when he first sets out, he doesn’t mind going back down the gradient to try finding a better spot. There’s a pretty low chance that he’d be able to find the highest point so early in the morning and he accepts that, contentedly moving back downhill before going higher. But as the day wears on, Dave gets a bit tired of the whole business. He’s passed a few peaks in the last few hours but none was as high as the one he passed around noon, and he’s wondering if maybe that midday mountain was the highest. So he becomes less happy about moving downhill and gradually starts just moving uphill. When night comes, he has found a good peak, but it’s not quite as high as the earlier one. He doesn’t much care about darkness, but he’s tired of searching, so he goes straight back to the highest point and makes camp, satisfied that this is most likely the best spot in the area.

In General

The general SA algorithm is analogous to Dave’s adventure. To apply SA to optimising an arbitrary system requires only two parts: a way to define the “energy” of the system’s state at any given time, and a temperature. We begin with a certain state s and a temperature T, and let’s say we’re looking for the state with minimum energy. We find a point (in configuration space) near s, called s'. The two states have energies e and e' respectively. If e' < e then s' is a better configuration than s, so we move to the new state and repeat from there. But if e' > e then we’d be moving in the wrong direction! This is where the temperature comes in: as the temperature of the simulation cools, the likelihood of moving to a “worse” state decreases. This is just like Dave getting tired through the day. So we define a function P(e,e',T) which compares the two energies, taking into account the temperature, and gives the probability of moving to the worse point. As the simulation evolves, the temperature cools, and the probability of choosing a worse configuration decreases.

The probability function P(e,e',T) can really be anything, as long as:

  • If e' < e (or > for a maximisation problem) then P = 1
  • If e' > e then P decreases as e'-e increases: there’s a lower chance of going up a steep gradient
  • P decreases as T decreases

In Particular

By keeping the truss’ member “network” (i.e. which members connect to which) constant, the state of the truss can be represented by the lengths of it’s members. The truss is represented by the (x,y) coordinates of it’s joints, so the list of joint positions acts as a proxy for the list of member lengths. Neighbouring states can then be found by moving a joint by a small amount; in this implementation I move a single joint by a random small amount each iteration. This way, the neighbouring states whose energies we compare differ only in the position of one joint.

The “energy” of a given truss is whatever we’re trying to optimise, which makes the maximum safe load a good metric. Another metric could be (maximum load)/(mass of truss) to give an efficiency score of sorts, but I haven’t explored this yet.

For the probability function I’ve found good results with P(e,e',T) = exp \left( -10 \cdot \frac{\Delta E}{T} \right). I’ve been running all these simulations on a five year old netbook, so there may be other functions which give better results but this one finds a good result in reasonable time (

For the temperature cycle, I’m repeatedly heating and cooling the simulation until sufficient time has passed since the last “best” state was found. The length of the cooling period is denoted by k_{max} – the maximum number of iterations k before reheating. The temperature function T(k,k_{max}) = \frac{2}{1+exp\left(15\cdot\frac{k}{k_{max}}\right)} has given decent results. It follows the lower half of a sigmoid curve, decreasing quickly from one and asymptotically approaching zero.

To limit the time taken, after 400 consecutive failed attempts at improving the truss, the simulation resets to the best solution so far. After 5 consecutive resets (2000 failed attempts in a row) the SA is finished. These numbers are arbitrary and as they’re quite small, generally won’t find the single optimum solution. But for a ten minute processing time, I find it sufficient.

As a side note, this program is an extension of the previous truss analysis script. In that script the truss was represented by a set of global tables, which was fine for that application but as SA requires comparing multiple trusses this wasn’t a great structure. So here the trusses are instances of a Truss class, which contains the necessary tables and methods for calculating the maximum load, etc. I’m more familiar with C++ classes so learning the python way was interesting!

And finally, here is another example of the program’s output:


The project is available on GitHub.

Something New

A sneak peak of things to come…

The Method of Joints: Part 2

See Part 1.

With the theory established, the python script was simple enough to write. Prior to this project I hadn’t used numpy before, but after seeing some examples online involving matrix operations, I figured numpy was the perfect tool.

The script is split into 9 major steps:

  1. Read data from the input files
  2. Construct the members table, filled with the coefficients l_{ij} and m_{ij}
  3. Construct the coefficient matrix C and the vector of applied forces P
  4. Solve the matrix equation!
  5. Find the maximum safe load for the truss
  6. Present the results

All the relevant data are kept in a few separate files, which are imported into some ndarrays:

  • containing the (x,y) coordinates of each numbered joint, in mm
  • a list of members with their two associated joints
  • the two supported joints and whether it is a pin or roller support
  • the forces applied to the truss, and their locations
  • the physical properties of the material we’re using

The first four use numpy.genfromtxt() to form their arrays, but is imported into a dictionary:

material = {}
with open(inpath+'', 'r') as fin:
    for line in fin:
        if (line[0] != '#'):
            spl = line.split()
            if len(spl)==2:
                (key, value) = spl 
                material[key] = float(value)

This way the material properties can be accessed easily, for example the Young’s modulus is material['E']

The members table initially has three columns: the member number, the number of joint i, and the number of joint j. The table is extended with a few empty columns, which are filled in with the relevant geometric properties of the member (taken from the (x,y) coordinates of it’s two joints): dx, dy, L, l_{ij}, m_{ij}.

Constructing C is pretty much as described in Part 1: the l_{ij} and m_{ij} values for each member are inserted in the appropriate place, and then the coefficients for the reaction forces are added in to the final three columns. These coefficients are \pm1 as they’re aligned with the axes. The applied forces are inserted into a vector “P
Solving the matrix equation is the easy part and it’s why I chose numpy in the first place:

Q = numpy.linalg.solve(C,P)

I’ve been setting up the applied forces to sum to 1, so Q now contains the portion of the load which is carried by each member. For example if Q[3] = -0.707, this means member #3 is carrying a compressive force of 0.707*(Load).

The maximum load for the truss as a whole is found by iterating through the members and finding the maximum load for the bridge before that member fails, then taking the lowest of these failure limits. This isn’t a rigorous analysis and only yielding and Euler buckling are accounted for, but that’s good enough for a bridge made of plastic straws:

maxL_t = np.zeros([members.shape[0],3])
maxL = 1e6
for k in range(members.shape[0]):
    if Q[k] != 0:
        # yield: F/A &lt;= (safety factor) * s_y
        maxL_t[k,0] = material['sf'] * material['s_y'] * material['A'] / abs(Q[k])
        # buckling: F &lt;= (safety factor) * pi^2 * E * I / l^2
        if Q[k] &lt; 0:
            maxL_t[k,1] = material['sf'] * numpy.pi**2 * material['E'] * material['I'] \
                            / (abs(Q[k]) * members[k,5]**2)
            # which is the smaller limit?
            maxL_t[k,2] = min(maxL_t[k,0], maxL_t[k,1])
            maxL_t[k,2] = maxL_t[k,0]
        if maxL_t[k,2] &lt; maxL:
            maxL = maxL_t[k,2]

And with that the analysis is complete. I’ve set up two modes of output: a plaintext mode, where the results are sent straight to a text file; and a presentation mode, where matplotlib is used to draw the bridge and the results are written to an .html page like this one!

The Method of Joints: Part 1

In one of my courses we’ve been tasked with building a truss bridge out of plastic tubes, slightly thicker than drinking straws, which must span 1 metre and hold 5kg of load. My team had a few basic designs that we wanted to start off with, which meant a whole bunch of calculation to see which one was best. Of course I’m not one to pass up an opportunity like that, and solving systems of linear equations is something that computers are kind-of alright at, so I wrote a python script that solves a given truss and returns the internal forces in each member. With the script complete, a new bridge design takes only about 2 minutes from data input to final results! The slowest step is sketching the bridge design…



The method of joints involves stepping through each joint in a truss and solving the equations of static equilibrium, i.e.

    \begin{align*} \sum{F_x} &= 0,\\ \sum{F_y} &= 0. \end{align*}

These two equations must be solved for each joint, giving 2n equations, where n is the number of joints present. The procedure is quite straightforward but having more than a few joints makes it get tedious very quickly. So, let’s make it faster.

A simple truss

Consider the frame above (with members 1 & 2 the same length). We have six unknowns to solve for: three reaction forces from the supports, and the three internal forces. We also have six equations describing the system, how convenient. These are:

    \begin{align*} 1:\sum{F_x} &= R_{1x} + F_3 \cdot \frac{1}{\sqrt{2}} = 0 \\ 1:\sum{F_y} &= R_{1y} - F_1 - F_3 \cdot \frac{1}{\sqrt{2}} = 0 \\ 2:\sum{F_x} &= R_{2x} + F_2 = 0 \\ 2:\sum{F_y} &= F_1 = 0 \\ 3:\sum{F_x} &= -F_2 - F_3 \cdot \frac{1}{\sqrt{2}} = 0 \\ 3:\sum{F_y} &= F_3 \cdot \frac{1}{\sqrt{2}} - 100N = 0 \\ &\Rightarrow F_3 \cdot \frac{1}{\sqrt{2}} = 100N \end{align*}

Now this system is simple enough to solve by inspection. Ignoring this though, we can represent this system as a matrix equation C \cdot F = P, where P is a vector of the applied (or known) forces, F is a vector of our unknowns, and C is the coefficient matrix. Thus:

    \begin{equation*} \begin{bmatrix} 0 & 0 & \frac{1}{\sqrt{2}} & 1 & 0 & 0 \\[0.3em] -1 & 0 & -\frac{1}{\sqrt{2}} & 0 & 1 & 0 \\[0.3em] 0 & 1 & 0 & 0 & 0 & 1 \\[0.3em] 1 & 0 & 0 & 0 & 0 & 0 \\[0.3em] 0 & -1 & -\frac{1}{\sqrt{2}} & 0 & 0 & 0 \\[0.3em] 0 & 0 & \frac{1}{\sqrt{2}} & 0 & 0 & 0 \end{bmatrix} \cdot \begin{bmatrix} F_1 \\[0.3em] F_2 \\[0.3em] F_3 \\[0.3em] R_{1x} \\[0.3em] R_{1y} \\[0.3em] R_{2x} \end{bmatrix} = \begin{bmatrix} 0 \\[0.3em] 0 \\[0.3em] 0 \\[0.3em] 0 \\[0.3em] 0 \\[0.3em] 100 \end{bmatrix} \end{equation*}

As long as the matrix C is non-singular, we can find it’s inverse and recover the solution F = C^{-1} \cdot P. In this case:

    \begin{equation*} \begin{bmatrix} F_1 \\[0.3em] F_2 \\[0.3em] F_3 \\[0.3em] R_{1x} \\[0.3em] R_{1y} \\[0.3em] R_{2x} \end{bmatrix} &= \begin{bmatrix} 0 & 0 & \frac{1}{\sqrt{2}} & 1 & 0 & 0 \\[0.3em] -1 & 0 & -\frac{1}{\sqrt{2}} & 0 & 1 & 0 \\[0.3em] 0 & 1 & 0 & 0 & 0 & 1 \\[0.3em] 1 & 0 & 0 & 0 & 0 & 0 \\[0.3em] 0 & -1 & -\frac{1}{\sqrt{2}} & 0 & 0 & 0 \\[0.3em] 0 & 0 & \frac{1}{\sqrt{2}} & 0 & 0 & 0 \end{bmatrix}^{-1} \cdot \begin{bmatrix} 0 \\[0.3em] 0 \\[0.3em] 0 \\[0.3em] 0 \\[0.3em] 0 \\[0.3em] 100 \end{bmatrix} = \begin{bmatrix} 0 \\[0.3em] -100 \\[0.3em] 141.4 \\[0.3em] -100 \\[0.3em] 100 \\[0.3em] 100 \end{bmatrix} \end{equation*}

That was easy! Now comes the question of somehow constructing the matrix for an arbitrary truss, without actually writing the equations by hand.

First, note that for a simple truss (constructed only of triangles) with n joints, there are 2n-3 members connecting them. If the truss is supported by three support forces, for example R_{1x}, R_{1y}, R_{2x} above, then the system is statically determinate: there are 2n equations describing the system and 2n unknowns to solve for.

An arbitrary member

Now consider an arbitrary member k in a truss, connected between two joints i and j. The equations

    \begin{equation*} i:\sum{F_x}, \quad i:\sum{F_y}, \quad j:\sum{F_x}, \quad j:\sum{F_y} \end{equation*}

are the only equations with any contribution from member k (as a member can only connect to two joints), and we can number these equations as 2i-1, 2i, 2j-1, and 2j respectively. The member k will appear in exactly these four equations, although the coefficient may be zero as with F_1 in the example above.

Looking at the x-component of forces at joint i, (equation number 2i-1):

    \begin{equation*} \dots + cos \theta \cdot F_k + \dots + P_{ix} = 0. \end{equation*}

We can simplify this a little by using ratios of the member’s dimensions rather than it’s angle. By defining

    \begin{align*} l_{ij} = cos \theta = \frac{\Delta x}{L_k} \\ m_{ij} = sin \theta = \frac{\Delta y}{L_k} \end{align*}

and by moving the known/applied force P_ix to the right hand side of the equation, we get the following:

    \begin{alignat*}{4} & \dots + l_{ij} \cdot F_k + \dots && = && -P_{ix} \quad && (eqn. \; 2i-1) \\ & \dots + m_{ij} \cdot F_k + \dots && = && -P_{iy} \quad && (eqn. \; 2i) \\ & \dots - l_{ij} \cdot F_k + \dots && = && -P_{jx} \quad && (eqn. \; 2j-1) \\ & \dots - m_{ij} \cdot F_k + \dots && = && -P_{jy} \quad && (eqn. \; 2j) \end{alignat*}

In the matrix C, the entire column k will be zero except for the rows corresponding to these equations. So the contributions to C from member k are:

    \begin{alignat*}{3} & C_{2i-1,k} \quad && = \quad && l_{ij} \\ & C_{2i,k} \quad && = \quad && m_{ij} \\ & C_{2j-1,k} \quad && = \quad && -l_{ij} \\ & C_{2j,k} \quad && = \quad && -m_{ij} \end{alignat*}

With C complete, we just construct the vector F as above, with F_k in order, and the three support reaction forces at the end. The vector of applied forces P is simple enough: for a load L applied at node i we just set

    \begin{alignat*}{3} & P_{2i-1} && = && -L_x \\ & P_{2i} && = && -L_y \end{alignat*}

And now we just find C^{-1}, solve F = C^{-1} \cdot P, and we’re done!