Over the course of several blog posts I will be documenting my progress in developing a very small series elastic actuator (SEA) for use in a robotic prosthetic hand.

Traditional actuators have a motor, gearbox and encoder making them very easy to use

]]>Over the course of several blog posts I will be documenting my progress in developing a very small series elastic actuator (SEA) for use in a robotic prosthetic hand.

Traditional actuators have a motor, gearbox and encoder making them very easy to use in robotics projects that require position or velocity control. However, while the output torque of DC motors can be controlled by monitoring their input current, it is quite difficult to do in the real world (as I discovered while working on an inverted pendulum). But by adding an elastic element (usually a spring) and encoder/potentiometer after the gearbox one can create a SEA that can measure and control the output torque (Hooke's law). Furthermore, SEA have the desirable property of compliance.

Over 3 years ago I designed and built my first robotic prosthetic hand, and to be frank, it was pretty terrible. But it was an excellent learning opportunity, amongst the many valuable lessons I discovered the importance of torque control and compliance in robotic systems. In the case of a robotic hand torque control is very important when picking up different objects; an egg requires a softer touch than a hammer. Additionally, compliance in a robotic hand is useful; if the hand was to impact something at a moderate force it is better to slightly and temporarily deform than to break on impact. Enter series elastic actuators, a brilliant type of actuator that resolves these two problems.

My unique design of a SEA is aimed to be small in size and cheap (<$60). The SEA will be designed to be about half the size of a whiteboard marker, such that it can fit inside the palm of the prosthetic hand, inline with its respectively finger (not including the thumb). The first image (above) shows a CAD model of the assembled SEA and below is an expanded view of the individual components.

From right to left the components in the image above are:

- Magnetic disk
- Hall effect encoder circuit
- N20 micro-gear motor (multiple gear ratios are available)
- Blue input frame: machined aluminium with axial hole for motor shaft, radial holes for set screw and torsion spring leg
- Black torsion spring with legs bent radially inwards
- Flange bearing with 3mm inner diameter for mounting on motor shaft and supporting the red output frame
- Red output frame: machined aluminium with axial holes for a 3x11mm rod and flange bearing, radial holes for set screw and torsion spring leg
- Flat disk potentiometer
- 3x11mm rod
- Bearing with 8mm inner diameter
- Bevel gear with set screw
- Output bevel gear with set screw

By measuring and calculating the angular difference between the magnetic encoder and potentiometer the spring force can be calculated $F_s=-k\theta$, thus enabling torque control when picking up objects. Note that the potentiometer does limit the output rotation to less that 180 degrees, however since this is designed to be used in a robotic hand this is not a concern. The output bevel gear will be the knuckle of the finger, directly rotating the proximal phalange and a linkage mechanism will move the distal and intermediate phalanges. If this SEA was to be adapted to a continuous rotation application I would trade the potentiometer for a gearbox and encoder (a gearbox will give more precise readings from the encoder).

This concludes the overview and scope for this project, I have begun manufacturing and am currently waiting on some supplies to arrive. There will be subsequent blog posts as I continue to develop this project.

]]>This tutorial will introduce a method of making a Christmas ornament that is dynamic and mesmerizing to watch. By spinning a rheoscopic fluid in a clear ornament ball you will be able to see the turbulent flow of the fluid inside. "Rheoscopic" means to ability to "see fluid currents"; this

]]>This tutorial will introduce a method of making a Christmas ornament that is dynamic and mesmerizing to watch. By spinning a rheoscopic fluid in a clear ornament ball you will be able to see the turbulent flow of the fluid inside. "Rheoscopic" means to ability to "see fluid currents"; this is fascinating optical phenomena is achieved by using a special fluid that contains light reflecting particles that do not dissolve in solution.

This fluid mixture is one of those simple yet captivating phenomena like fire and magnetism that entertains children and adults alike.

- Water (tap water works, distilled water is best)
- Rheoscopic concentrate (I used Pearl Swirl, but there are multiple options)
- Food colouring
- Clear, fillable Christmas ornaments (can be found at dollar stores and craft stores)
- Hot glue
- Coreless DC motor with drone rotor (like this)
- Electrical wire, heat shrink tube, switch, fuse, solder, 3V power supply

In a container mix 1L of water, 1 tablespoon of rheoscopic concentrate and a couple drops of food colouring. As you begin to mix the solution you will see chaotic flows arise. For the best results I recommend boiling the water first (to remove dissolved air), I personally didn't do this but I think the final result would look better if you did.

Note that if the fluid is allowed to rest for several hours most of the light reflecting particles will settle at the bottom of the container. This isn't an issue as long as you scrap the very bottom of the container when you come back to stir it. I left my own mixture for several days and within 20 seconds of mixing the rheoscopic effect came back fully.

At this point it is clear that a stir rod works very well for stirring the mixture. But I encourage you to experiment with different methods as it is a lot of fun to watch the turbulent flows. Some ideas that were successful for me: dropping objects in the container, a magnetic stir plate (as found in chemistry classrooms/labs), a tube to blow air through, and placing the room temperature fluid on a hot surface (to visualize thermal convection currents).

But coming back to the objective of making Christmas ornaments with this fluid, I found the best method was to use a small DC motor to stir the mixture (after all who want to aggressively rattle ornaments)?

It may seem silly at first to put a DC motor inside a fluid, and while generally it is not a great idea, for short term applications (like the couple of weeks until Christmas) it will not lead to an immediate failure. As such we can submerge the motor directly in the fluid to agitate it. Nonetheless, the likely cause of failure will still be corrosion from electrolysis on the brushes. I will come back and update this when the motors do fail, so far they have been going for 48 hours with no problems. Using distilled water will decrease the chance of failure compared to regular tap water.

I came up with 2 different configurations:

The vibration motor is commonly used in mobile phones as a notification feature. By soldering thin wires are submerging the motor it can agitate the rheoscopic fluid, although the amplitude of agitation is small. The better performing motor was the coreless DC motor with drone rotor blade attached.

Furthermore, the coreless DC motors only has open holes at their base (which will be sealed with hot glue), the end of the motor that has the axle coming out tightly mates with a metal bushing, making it difficult for any water to ingress, thus further delaying the potential electrical failure.

When the silver top of the clear, empty ornament is pulled off there will be two holes that the wires from the motor can be fed through (see below pictures). I then tied a reef knot to temporarily secure the motor in place. To permanently secure the motor I applied hot glue to the base of the motor, do your best to fully cover the base of the motor as this is the most likely point that fluid will get inside the motor.

Before attaching the rotor blades you will need to bend them to fit through the neck of the ornament. I was able to bend mine without issue but in case yours is made of a brittle plastic you may want to heat them with a hot air gun first. Once bent, attach them to the motor (see pictures).

To transfer the fluid to the ornaments I found it easiest to use a plastic syringe but a funnel would work too. Remember to stir the mixture before filling the ornaments to ensure none of the particles have settled at the bottom of the container. Leave a little bit of air in the neck of the ornament as the motor will displace some of the fluid.

Fill the silver cap of the ornament around the motor with hot glue and quickly and firmly push it onto the neck of the ornament, ensuring that a waterproof seal is made between the surfaces.

Before proceeding to the next step ensure each ornament works by connecting the motor terminals to a 3V power source.

Now that the ornaments have been made there are a few final steps to make them ready to hang. Firstly cut 2 lengths of wire with openings every 30 centimetres or so, at each interval solder a connection to to the motor. I personally used some male-female connections (link) so that if one of the motors fails I can easily remove it without unsoldering everything. Add heat shrink tubing to insulate the solder joints, preventing short circuits. I also recommend connecting a fuse so that should one ornament fail the power supply will disconnect instead of melting the wires/damaging the power supply.

Finally I wrapped the wires connecting the ornaments with some tinsel (you could also use LED lights if you wish). Hook the wires up to a 3V power source and you have a chain of ornaments, ready to hang.

Suggested improvement:

- One improvement that I may add in the next version is instead of a fixed power supply I might try using a PWM driver to control the speed and direction of the motor. Personally I like the slower speeds for visualizing the fluid currents.

- Derivation of Simple Pendulum (Python Simulation) - link

- Derivation of Simple Pendulum (Python Simulation) - link
- Building a Physical Inverted Pendulum
- Controlling Inverted Pendulum with PID controller
- State Space Model of Inverted Pendulum - link
- Designing the Optimal Controller for an Inverted Pendulum (LQR control) - link
- Simulating an Inverted Pendulum with an LQR Controller (Python) - link
- Deploying an LQR Controller on a Physical Inverted Pendulum

The associated code can be found on GitHub and a video series on YouTube (not currently released as of Nov. 2020). While this project aims to be encompassing and self contained, it does assume familiarity with classical mechanics, linear algebra, differential equations, some control theory, and programming (Python).

In the previous 2 blog posts I have explored the state space model of an inverted pendulum and derived the optimal LQR controller for stabilizing the system. Before going to the trouble of building a physical inverted pendulum lets first verify that the controller works as expected in simulation.

If you wish to dive right into the code it can be found on my GitHub

To provide an explanation of the code line by line is difficult outside of video format. As such I will touch on the highlights. There are two files: variables.py and simulation.py. If you have read the previous 2 blog posts the variables.py file should be fairly clear, I have just defined all the relevant variables that have been previously discussed. Feel free to adjust the mass of the cart, rod, damping etc. You can also experiment with different values in the Q and R matrices, as discussed in the previous blog post.

The file simulation.py builds off the existing code from part one of this series, the animate() function is largely unchanged. The primary difference is that the non-linear equation for the simple pendulum has been replaced by the function StateSpace(t, i). This function evaluates the non-linear inverted pendulum with the applied LQR input control.

]]>- Derivation of Simple Pendulum (Python Simulation) - link

- Derivation of Simple Pendulum (Python Simulation) - link
- Building a Physical Inverted Pendulum
- Controlling Inverted Pendulum with PID controller
- State Space Model of Inverted Pendulum - link
- Designing the Optimal Controller for an Inverted Pendulum (LQR control) - link
- Simulating an Inverted Pendulum with an LQR Controller (Python) - link
- Deploying an LQR Controller on a Physical Inverted Pendulum

The associated code can be found on GitHub and a video series on YouTube (not currently released as of Nov. 2020). While this project aims to be encompassing and self contained, it does assume familiarity with classical mechanics, linear algebra, differential equations, some control theory, and programming (Python).

This blog post builds directly from my previous blog post on the state space model of an inverted pendulum. In this post I will derive the optimal controller for the inverted system using control theory, while I aim to cover all the relevant steps this is not designed to be a text book as such I will be presenting equations and theorems without mathematical proofs. For good online control theory courses I would recommend MIT OpenCourseWare and Steve Brunton's YouTube series on control theory.

To recap from the previous blog post we found that the state space model is:

$$\begin{bmatrix} \dot{x} \\ \dot{\theta} \\ \ddot{x} \\ \ddot{\theta} \end{bmatrix} = \begin{bmatrix} 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & \frac{m_2^2gl^2}{m_1m_2l^2 + I(m_1 + m_2)} & \frac{- \beta(m_2 l^2+I)}{m_1m_2l^2 + I(m_1 + m_2)} & 0 \\ 0 & \frac{m_2gl(m_1+m_2)}{m_1m_2l^2 + I(m_1 + m_2)} & \frac{-m_2l \beta}{m_1m_2l^2 + I(m_1 + m_2)} & 0 \end{bmatrix} \begin{bmatrix} x \\ \theta \\ \dot{x} \\ \dot{\theta} \end{bmatrix} + \begin{bmatrix} 0\\0\\ \frac{m_2l^2+I}{m_1m_2l^2 + I(m_1 + m_2)} \\ \frac{m_2l}{m_1m_2l^2 + I(m_1 + m_2)} \end{bmatrix} u$$

$$y=\begin{bmatrix} 1 & 0&0&0\\0&1&0&0\end{bmatrix}\begin{bmatrix} x \\ \theta \\ \dot{x} \\ \dot{\theta} \end{bmatrix} +\begin{bmatrix} 0 \\ 0\end{bmatrix}u$$

Which, more succinctly, can be expressed in the general form:

$$\dot{\vec{x}} = A\vec{x}+Bu$$

$$\vec{y} = C\vec{x}+Du$$

To ensure our efforts are not futile it is a good first step to ensure the system we wish to control is controllable to begin with. A system is controllable if for all $\vec{x}_0$, $\vec{x}_1$, and $T>0$ there exists a control input, $u$, such that $\vec{x}(t=0)=\vec{x}_0$ and $\vec{x}(t=T)=\vec{x}_1$. The test for controllability states that if the rank of the matrix $\begin{bmatrix} B & AB & ... & A^{n-1}B\end{bmatrix}$ is $n$ (where A is an $n \times m$ matrix ) then the system is controllable (proof omitted).

The following is controllability matrix evaluated for the following quantities: $m_1=1kg, m_2=0.25kg, g=9.81m \cdot s^{-2}, l=0.5m, \beta = 1 kg \cdot s^{-1}$.

$$\begin{bmatrix} B & AB & ... & A^{n-1}B\end{bmatrix} = \begin{bmatrix} 0 & 0.94 & -0.89 & 3.28 \\ 0 & 1.41 & -1.33 & 25.69 \\ 0.94 & -0.89 & 3.28 & -5.39 \\ 1.41 & -1.33 & 25.69 & -27.63\end{bmatrix}$$

I used the Python command *control.ctrb(A, B)* to calculate this matrix and the numpy.linalg.matrix_rank command to confirm that the rank of the matrix is 4. Thus the system is controllable.

While we can intuitively tell that the inverted pendulum is an unstable system, we should still test for it mathematically, in some systems that are not as trivial as the inverted pendulum it is not immediately obvious from inspection if the system is stable or not. There are several tests for stability but I will use the eigenvalue method here. The proof is omitted from this blog post but a system is internally stable if the all the real parts of the matrix A are negative. That is to say the system is stable iff $Re(\lambda)<0, \quad \forall \lambda$.

Using the Python command, scipy.linalg.eig(A), the eigenvalues are found to be:

$$\begin{bmatrix} 0+0i \\ 4.1 + 0i \\ -4.2+0i \\ -0.8 +0i \end{bmatrix}$$

Thus since some of the eigenvalues have non-negative real parts the system is unstable.

Now that it has been established that the system is controllable but unstable, we can find a suitable control input, $u$, to stabilize the system. Lets define the control input in terms of the state vector, $\vec{x}$:

$$u=-K(\vec{x}(t) - \vec{x}_{setpoint})$$

Where $\vec{x}_{setpoint} = \begin{bmatrix}0& \pi & 0 &0 \end{bmatrix}^T$ is our desired final state and $K$ is some yet to be defined matrix. Notice that our state space can be rewritten as follows:

$$\dot{\vec{x}} = A\vec{x}+Bu = (A-BK)\vec{x}+K\vec{x}_{setpoint}$$

We now have the freedom to choose the matrix K such that the eigenvalues of $(A-BK)$ lead to stability, as shown in the previous section. Let me choose (somewhat) randomly the matrix values $K = \begin{bmatrix} -100 & 800 & -150 & 200\end{bmatrix}$. If we repeat the stability test for $(A-BK)$ we get the following eigenvalues:

$$\begin{bmatrix} -134.7+0i \\ -4.6 + 0i \\ -1.4+0.5i \\ -1.4 -0.5i \end{bmatrix}$$

Notice that now all the eigenvalues have negative real parts, thus the system is stable. This is great! We now have a control input, $u = -K(x-\vec{x}_{setpoint})$, that if applied to the inverted pendulum system would stabilized. However we can do better, we want to find the optimal control matrix K for the most efficient control.

The procedure to do this is a little involved and to be honest I don't want to write an entire control theory book. As such I will present the relevant equations and an overview of the steps required to find the optimal matrix K. If you want a deeper understanding of what these equations mean, how they work, and the associated mathematical proofs I advise consulting online control theory courses like the one available through MIT OpenCourseWare.

We can define a cost function $J$ as follows. It will be our task to minimize this cost function, in doing so we will find the optimal control.

$$J(\vec{x}, u) = \int_{0}^{\infty}(x^TQx+u^TRu)dt$$

Where $Q$ is a positive semi-definite matrix and $R$ is positive-definite, both are symmetric. $Q$ is effectively used to express the cost we wish to place on physical system: i.e. how important is it to get the correct final position, the speed and aggression at which to move? $R$ is effectively the cost we wish to apply to the input, in this case the motor, maybe we wish to conserve power or perhaps we have a powerful motor and don't care how hard we force it to move. A common starting point is to define them as:

$$Q=C^TC=\begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 &1 & 0 \\ 0& 0& 0& 1 \end{bmatrix} \quad R =\begin{bmatrix} 1 \end{bmatrix} $$

If we wanted to prioritize the final angle, $\theta$, and had a powerful motor to use we could adjust these to:

$$Q=\begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 10 & 0 & 0 \\ 0 & 0 &1 & 0 \\ 0& 0& 0& 1 \end{bmatrix} \quad R =\begin{bmatrix} 0.001 \end{bmatrix} $$

Typically these variables are changed by factors of 10. In the next blog post (simulating an inverted pendulum in Python) we will be able to experiment with various $Q$ and $R$ matrices to see the different effects they can have.

By the minimization of the quadratic cost function, $J(\vec{x}, u)$, the optimal controller, $K$, can be expressed as:

$$K = R^{-1}B ^TP$$

Where P is a symmetric positive-definite matrix found from the algebraic Ricatti equation (derivation is omitted):

$$A^TP+PA-PBR ^{-1}B^TP +Q =0$$

To calculate the matrix K this in Python the following commands can be used:

P = np.matrix(linalg.solve_continuous_are(A, B, Q, R))

K = np.matrix(linalg.inv(R)*(B.T*P))

Thus we have found the optimal control input:

$$u = -K(\vec{x}-\vec{x}_{setpoint}) = -R^{-1}B ^TP(\vec{x}-\vec{x}_{setpoint})$$

Finally we can once again check the eigenvalues of $(A-BK)$ and confirm that their real parts are negative:

$$\begin{bmatrix} -53.8+0i \\ -1.1 + 0i \\ -2.8+0.2i \\ -2.8 -0.2i \end{bmatrix}$$

For the last two blog posts I have been referring to the magical mathematical entity $u$, but what is it in the real world? Recall in the second section of the previous blog post I chose to redefine the force applied to the cart, $\vec{F}$, as the control input, $u$. This was predominately to adhere to the standard control theory conventions, if we now go back to the force, we can express the force applied to the cart as the torque, $\tau$, from the motor and gearbox:

$$F=\frac{\tau}{Gr} \quad where \enspace G=\frac{N_2}{N_1}$$

Where $r$ is the armature radius, $G$ is the gear ratio, and $N_1, N_2$ are the number of teeth on the first and second gears, respectively (assuming a 2 piece gear train).

Making the substitution the control input can be expressed in terms of the motor torque, something that we can control via hardware and software by regulating the current applied to the motor.

]]>The cube is made of 6 identical* face pieces that are designed to interlock with each other; each face has 9 acrylic tiles for diffusing the LED light and 18 brass square rods that act as capacitive touch sensors. Behind each of these faces is a circuit board with 9 different RGB addressable LEDs. The circuit design is discussed further in the next section.

The acrylic tiles are 3.175mm thick with a side length of 15.5mm; they are supported in the 3D printed frame through a friction press and a small 0.5mm lip in the 3D printed frame to ensure they are not pushed in too deep. The edge acrylic pieces have their sides tapered at 45 degrees (see pictures), the current design utilizes friction to hold them in place but they could be glued in place too.

The brass square rods are 3.175 x 3.175 x 15.5mm and also held in place through friction. Small channels are cut in the 3D printed frame for connecting the brass rods to the PCB via wires.

Inside the cube, the six cube-face circuit boards are connected to a seventh control board and communicate via an I2C bus. Behind the circuit boards is a smaller cube-shaped void, in this void a lithium-ion battery is placed.

In the first image note the small channels that connect the PCB board behind the frame to the brass capacitive touch sensors at the front of the frame. In the second image you can see a the 3D printed model with 9 acrylic tiles and 1 of the 18 brass capacitive touch sensors.

*modification is made to one of the faces to allow for a battery charging port

See block circuit diagram in the pictures for a broad overview of the circuit

Each of the 3D printed frames has a two-sided circuit board mounted on it. On one side will be 9 RGB addressable LEDs for displaying the colour of the Rubik's cube tiles. The 18 brass capacitive touch sensors will be connected to the circuit board via short wires that fit into the channels cut out of the 3D printed frame, these 18 signal lines will be split into two groups running to one of two MPR121 chips. The MPR121 is a 12 channel capacitive touch chip with four I2C addresses (datasheet link). Due to the limited number of I2C addresses and the need for 12 MPR121 chips (2 for each of the 6 cube faces) a central I2C switch is included for cycling through each of the 6 circuit boards. This central I2C switch is mounted on a seventh control board, that is positioned in the central void of the cube, next to the battery. In addition to controlling the I2C communication it will also have a micro-controller, voltage regulator, and charge control circuits.

The micro-controller will be responsible for periodically collecting the binary state of the 108 capacitive touch sensors (touched/not touched). If a row/column of capacitive touch sensors on any of the faces are touch one after the other and within a certain time threshold a “swipe” will be registered and the LEDs will change colour to reflect a ‘rotation’ of the cube. Touches that are detected that are not part of a “swipe” will be ignored since in order to hold the cube the user will inevitably be making contact with some of the capacitive touch sensors.

There are two noteworthy shortcomings in the current proof-of-concept design: cube size and battery life. The current design has a side length of 65.3mm whereas the standard Rubik's cube has a side length of 57mm. In hindsight, I could have made a design that adheres to this measurement but I have already cut all the acrylic pieces by hand (most of which are tapered at 45 degrees on one or two of their faces), and I would rather not go through this laborious task again. After all this is merely a fun, proof-of-concept project.

Secondly, the short battery life (currently) prevents this project from being little more than a novelty. I am using a 3.7V 1800mAh lithium-ion battery; this was the largest single cell battery I could find that would fit in the interior cavity of the cube. This is a concern as 54 RGB LEDs will consume significant power, rough estimates put the expected operation time at about 30 minutes. To increase the battery life in future versions I will explore the possibility of using multiple smaller battery cells wired in parallel to form a "battery cube", this design would optimize the efficient use of the interior space. I have not done this in my current design as wiring (and charging) lithium-ion batteries in parallel is non-trivial can be dangerous if not done correctly; I do not currently have the knowledge or experience to feel safe doing this.

Future developments could include installing an inertial measurement unit (e.g. shake the cube in frustration to reset to the solved state), a Bluetooth/phone connection (for tracking speeds, logging progress data, and competing against friends), and a learning program: the cube makes a move and waits for the user to replicate it, repeating this process until the cube is solved (used to teach the user different Rubik’s cube solving algorithms).

This build log documents the PCB design of the 6 identical boards for each of the 6 cube faces. Recall from the 'Project Details' section (above) that there will be a 7th central PCB that has a microcontroller and is responsible for power regulation. Each of the 6 cube-face boards have the following features:

- 9 WS2812B LEDs (datasheet)
- 2 MPR121 chips (datasheet)
- 18 through-holes for connecting brass capacitive touch sensors
- 4 mounting holes (M2 mounting bolts)
- Additional through-holes for inter-PCB cabling

Below are images of the top and bottom PCB traces (red and blue lines, respectively). Note that some traces appear to be touching, this is not the case; I merely used a screenshot and then poorly compressed the image.

The topside has the 9 LEDs chained together and can then be connected to the neighbouring board, once connected there will be a continuous strip of 54 addressable LEDs instead of 6 segments of 9 LEDs. Matrix rotations will be used for calculations in the software.

The bottom side has two 12-channel MPR121 chips with different I2C addresses of 0x5B and 0x5C. These chips were chosen as the I2C protocol allows for easily checking the state of individual capacitive touch sensors as opposed to a convoluted multiplexing network. This will make developing the software that monitors for "swipes" easier (and faster). I unfortunately could not source an alternative 18-channel capacitive touch chip with I2C. If you wish to replicate this project please be aware that this chip is a legacy device and not recommended for new designs, however since this is a hobby project and not something I plan to make commercially this was not a concern (available for purchase here-as of Oct. 2020). Furthermore, on a personal level, I wanted to practice reflow soldering with a QFN chip (quad-flat no-leads).

Note that 6 of the 18 through-holes for connecting the brass capacitive touch sensors are on the edge of the board as castellated holes. While not my first preference this was the only way I could devise that connected the PCB to the corner capacitive touch sensors without interfering or block the light from the LEDs.

The PCBs were manufactured by OSH Park.

Assembling the PCBs was done with reflow soldering using a hot air gun. I would recommend doing the bottom side first to get the QFN packages out of the way. Below is an image of the assembled board mounted on a test rig.

Please stand by. Currently optimizing PCB design based on initial testing.

]]>- Derivation of Simple Pendulum (Python Simulation) - link

- Derivation of Simple Pendulum (Python Simulation) - link
- Building a Physical Inverted Pendulum
- Controlling Inverted Pendulum with PID controller
- State Space Model of Inverted Pendulum - link
- Designing the Optimal Controller for an Inverted Pendulum (LQR control) - link
- Simulating an Inverted Pendulum with an LQR Controller (Python) - link
- Deploying an LQR Controller on a Physical Inverted Pendulum

The associated code can be found on GitHub and a video series on YouTube (not currently released as of Nov. 2020). While this project aims to be encompassing and self contained, it does assume familiarity with classical mechanics, linear algebra, differential equations, some control theory, and programming (Python).

A sketch of the inverted pendulum system is shown in the figure below. From this we will derive the state space model (and eventually find the optimal controller to balance the inverted pendulum). The overall process for finding the state space model will be as follows: find the Lagrangian, evaluate the Euler-Lagrange equations, rearrange the coupled differential equations in state space form, and then linearize by evaluating the system at $\theta=\pi$ (pendulum vertical position).

The Lagrangian is defined as $\mathcal{L}=T-V$, where $T$ and $V$ are the kinetic and potential energies, respectively. The terms in the sketch, above, are largely self evident, however I will choose to define $l=L/2$. I also make the assumption that the rod is made of a uniform mass, thus making $l$ the centre of mass of the rod. Using this information we can establish:

$$\dot{\vec{v}}_{cart} = \dot{x}\hat{i} $$

$$\dot{\vec{v}}_{rod} = \frac{d}{dt}\Big( x+l\sin(\theta) \Big)\hat{i} + \frac{d}{dt}\Big( l\cos(\theta) \Big)\hat{j} = (\dot{x} + l\dot{\theta}\cos(\theta))\hat{i} + (l\dot{\theta}\sin(\theta))(-\hat{j}) $$

$$T = \frac{1}{2}m_1v_{cart}^2 + \frac{1}{2}m_2v_{rod}^2 + \frac{1}{2}I_{rod}\omega^2 $$

$$V = m_2gl(1-\cos(\theta)$$

where $I_{rod}$ and $\omega$ are the moment of inertia and angular velocity of the rod, respectively. Substituting the above four equations into the Lagrangian, $\mathcal{L}=T-V$.

$$\mathcal{L}=\frac{1}{2}m_1\dot{x}^2 + \frac{1}{2}m_2(\dot{x}^2+\dot{x}l\dot{\theta}\cos(\theta) + l^2\dot{\theta}^2) + \frac{1}{2}I_{rod}\dot{\theta}^2-m_2gl(1-\cos(\theta)$$

From the Lagrangian by using the Euler-Lagrange equations, below, the equations of motion of the system can be found.

$$\frac{d}{dt}\Big( \frac{\partial\mathcal{L}}{\partial \dot{q_i}} \Big) - \frac{\partial\mathcal{L}}{\partial q_i} = Q_i \quad \quad \quad i =x,\theta$$

Where $Q_{\theta}, Q_x$ are the external forces applied to the system in the $\theta$ and $x$ axes. In the $x$-axis there is a force applied to the cart (that we will later use to stabilize the inverted pendulum) and some dampening force proportional to the velocity. As such $Q_x = F-\beta \dot{x}$. In the $\theta$-axis there is also some friction at the pivot point however this resistance is negligible. As such I will define: $Q_{\theta}=0$. These quantities can be substituted into the Euler Lagrange equations.

For the sake of brevity (and my own sanity), I'm opting out of transcribing the partial derivatives in LaTeX. I will leave it as an exercise to verify that when the Euler-Lagrange equations are evaluated for $i =x,\theta$ that the equations of motion come out to be:

$$m_2\ddot{x}l\cos(\theta)+(m_2l^2-I)\ddot{\theta}= -m_2gl\sin(\theta)$$

$$(m_1+m_2)\ddot{x} + m_2l\ddot{\theta}\cos(\theta) = m_2l\dot{\theta}^2\sin(\theta) + F - \beta\dot{x}$$

Going forward I will rename the force, $F$, as $u$ which is the convention in control theory for a 'control input'. In subsequent blog posts I will clarify this convention and explain how $F$ can be expressed in terms of the motor torque and how in the real world the torque is controlled by the electrical current to the motor. But for the moment we will live in mathematical bliss and leave the control input as simply $u$.

At this point we have found the coupled equations of motion of the system. However $\ddot{x}$ and $\ddot{\theta}$ appear in both equations. In this section I will rearrange these equations in matrix form and through some tedious manipulations separate the $\ddot{x}$ and $\ddot{\theta}$ terms.

In matrix form, the coupled equations of motion can be expressed in matrix form as:

$$\begin{bmatrix} m_2l\cos(\theta) & m_2l^2 +I \\ m_1+m_2 & m_2l\cos(\theta)\end{bmatrix} \begin{bmatrix} \ddot{x} \\ \ddot{\theta} \end{bmatrix} = \begin{bmatrix} 0 & 0\\ -\beta & m_2l\dot{\theta}\sin(\theta) \end{bmatrix} \begin{bmatrix} \dot{x}\\ \dot{\theta} \end{bmatrix} + \begin{bmatrix} 0 & 0 \\ 0 & 0 \end{bmatrix} \begin{bmatrix} x \\ \theta \end{bmatrix} + \begin{bmatrix} 0 & -m_2gl\sin(\theta) \\ 0 & 0 \end{bmatrix} \begin{bmatrix} 1 \\1 \end{bmatrix} +\begin{bmatrix} 0 \\ 1 \end{bmatrix} u$$

Refactoring the terms a more concise form is:

$$E\ddot{\vec{x}} = F\dot{\vec{x}} + G \vec{x} + H + Iu$$

And by taking the inverse of $E$:

$$\ddot{\vec{x}} = E^{-1}F\dot{\vec{x}} + E^{-1}G \vec{x} + E^{-1}H + E^{-1}Iu$$

$$E^{-1} = \frac{1}{m_2^2l^2 \cos^2(\theta) - (m_1 + m_2)(m_2l^2+I) } \begin{bmatrix} m_2l\cos(\theta) & -m_2l^2 -I \\ -m_1-m_2 & m_2l\cos(\theta)\end{bmatrix}$$

By doing some tedious matrix multiplication (exercise for the reader) the expression for vector $\ddot{\vec{x}}$ can be expanded to:

$$\ddot{x} = \frac{1}{D} \Big[ \beta(m_2l^2 +I)\dot{x} - m_2l\dot{\theta}^2\sin(\theta)(m_2l^2 +I) - m_2^2gl^2\cos(\theta)\sin(\theta) - (m_2l^2 +I)u\Big]$$

$$\ddot{\theta} = \frac{1}{D} \Big[ -\beta m_2l\cos(\theta)\dot{x} + m_2^2l^2\dot{\theta}^2\cos(\theta)\sin(\theta) +(m_1 + m_2)m_2gl\sin(\theta) + m_2l\cos(\theta)u \Big]$$

$$D = m_2^2l^2 \cos^2(\theta) - (m_1 + m_2)(m_2l^2+I)$$

To recap, so far we have found the equations of motion of the system and expressed them as separate equations of $\ddot{x}$ and $\ddot{\theta}$. These non-linear equations will be useful in the blog post regarding simulating the physical system. However, to find the optimal LQR controller these equations must be linearized around the unstable equilibrium point $\theta = \pi$. Since the LQR controller is only valid for near the equilibrium point we can assume that the system remains within $\phi$ of $\pi$, where $\phi$ is a small angle. Therefore the nonlinear terms of interest become:

$$\cos(\theta) = \cos(\pi + \phi) = \cos(\pi) - \sin(\pi\phi) + O^2 \approx -1$$

$$\sin(\theta) = \sin(\pi + \phi) = \sin(\pi) + \cos(\pi\phi) + O^2 \approx -\phi$$

$$\dot{\theta}^2 = \dot{\phi}^2 \approx 0$$

And so, using this information $\ddot{x}$ and $\ddot{\theta}$ become:

$$\ddot{x} = \frac{1}{D} \Big[ \beta(m_2l^2 +I)\dot{x} - m_2^2gl^2\theta - (m_2l^2 +I)u\Big]$$

$$\ddot{\theta} = \frac{1}{D} \Big[ \beta m_2l\dot{x} -(m_1 + m_2)m_2gl\theta - m_2lu \Big]$$

$$D = -m_1m_2^2l^2 - (m_1 + m_2)I$$

To review, the state space model representation is a system of two equations (below). Where A, B, C, and D are matrices that describe the state, input, output, and feedthrough, respectively. In most cases, including the inverted pendulum, D is a zero matrix.

$$\dot{\vec{x}} = A\vec{x}+Bu$$

$$\vec{y} = C\vec{x}+Du$$

Looking at the pair of equations at the end of part 4, we can see that by converting the linear system back into matrix form this will yield the general state space form:

$$\begin{bmatrix} \dot{x} \\ \dot{\theta} \\ \ddot{x} \\ \ddot{\theta} \end{bmatrix} = \begin{bmatrix} 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & \frac{m_2^2gl^2}{m_1m_2l^2 + I(m_1 + m_2)} & \frac{- \beta(m_2 l^2+I)}{m_1m_2l^2 + I(m_1 + m_2)} & 0 \\ 0 & \frac{m_2gl(m_1+m_2)}{m_1m_2l^2 + I(m_1 + m_2)} & \frac{-m_2l \beta}{m_1m_2l^2 + I(m_1 + m_2)} & 0 \end{bmatrix} \begin{bmatrix} x \\ \theta \\ \dot{x} \\ \dot{\theta} \end{bmatrix} + \begin{bmatrix} 0\\0\\ \frac{m_2l^2+I}{m_1m_2l^2 + I(m_1 + m_2)} \\ \frac{m_2l}{m_1m_2l^2 + I(m_1 + m_2)} \end{bmatrix} u$$

$$y=\begin{bmatrix} 1 & 0&0&0\\0&1&0&0\end{bmatrix}\begin{bmatrix} x \\ \theta \\ \dot{x} \\ \dot{\theta} \end{bmatrix} +\begin{bmatrix} 0 \\ 0\end{bmatrix}u$$

]]>- Derivation of Simple Pendulum (Python Simulation) - link

- Derivation of Simple Pendulum (Python Simulation) - link
- Building a Physical Inverted Pendulum
- Controlling Inverted Pendulum with PID controller
- State Space Model of Inverted Pendulum - link
- Designing the Optimal Controller for an Inverted Pendulum (LQR control) - link
- Simulating an Inverted Pendulum with an LQR Controller (Python) - link
- Deploying an LQR Controller on a Physical Inverted Pendulum

This blog is supplemental documentation pertaining to the derivation for the equation of motion used in the python program, simple-pendulum.py

While the program uses the scipy.integrate.solve_ivp() function, I will also touch on the small angle approximation for solving the system analytically for small initial angles $\theta_0 << 1$. There are a few different ways to derive the equations of motion, I will be using the torque derivation, and mention in passing the Lagrangian approach.

Firstly, since this is a rod-pendulum with uniform mass (unlike the slightly more common, point mass attached to a massless rod pendulum), the moment of inertia $I$ must be found.

$$I_{rod,CoM} = \int_{0}^M r^2dm = \int_{-\frac{L}{2}}^{\frac{L}{2}}\frac{r^2m}{L}dr = \frac{1}{12}mL^2 $$

Using the parallel axis theorem, we can shift the moment of inertia from the center of mass (CoM) to one of the ends:

$$I_{rod,end} = I_{rod,CoM} + md^2 = \frac{1}{3}mL^2$$

We now have the information required to find the equation of motion via the torque:

$$ \tau = I\alpha $$

$$ -mg\sin(\theta)(\frac{L}{2}) = I_{r,e}\frac{d^2 \theta}{d t^2}$$

$$\frac{d^2 \theta}{d t^2} = - \omega^2 \sin(\theta)$$

where $\omega = \sqrt{\frac{3g}{2L}}$ is the natural frequency of the pendulum. At this point if using a small initial angle $\theta_0 << 1$, one could solve the system analytically using the small angle approximation $\sin(\theta) \approx \theta$, and then integrating twice to give $\theta(t) = \theta_0 \cos(\omega t)$. Refer to the first couple of commits in the repository to see the code for a small angle approx (of a point-mass pendulum). At this point the pendulum animation can be created using matplotlib. Note: this blog post focuses on the physics of the system, not the actual software. For a deeper understanding of how scipy.integrate.solve_ivp() is used please refer to the scipy documentation.

To create the phase space plot I will calculate the Hamiltonian and Lagrangian (optional) of the system. The kinetic and potential energy can be found as:

$$T = \frac{1}{2}mv^2 + \frac{1}{2}I\dot\theta^2 = \frac{1}{6}mL^2\dot\theta^2$$

$$V = \frac{mgL}{2}(1-\cos(\theta))$$

The Lagrangian, $L=T-V$, can be used with the Euler-Lagrange equation to find the equation of motion (as opposed to the torque derivation, above). This is left as an exercise for the reader. Using the generalised coordinate, $q=\theta$, and generalized momentum, $p = \frac{\partial L}{\partial \dot\theta}$, the Hamiltonian can be expressed as:

$$H = T + V = \frac{p^2}{6m L^2} + mgL(1-\cos(q))$$

The Hamiltonian represents the total energy of the system and the different energy level contours are expressed in the phase space plot.

]]>Using a 1kg load cell, HX711 amplifier chip and an ATmega328P I created the following battery operated kitchen scale. While some DIY projects help save money, this is not one of them. Cheap kitchen scales can be found for approximately $15 (CAD) and the individual parts for this project total about $50 (CAD). Nonetheless for those tinkerers who enjoy making stuff for the sake of making it the following video tutorial may be useful along with the wiring diagram and programs, down below. Link to GitHub repo.

```
#include "HX711.h"
HX711 scale;
#define CLK A0
#define DOUT A1
int calibration_factor = 1000;
void setup() {
Serial.begin(9600);
Serial.println("Kitchen scale calibration");
Serial.println("Remove all weight from scale");
scale.begin(DOUT, CLK);
scale.set_scale();
scale.tare(); //Reset the scale to 0
Serial.println("Place an object of known weight on the load cell");
Serial.println("Adjust calibration factor using the following keys");
Serial.println("'a' or 'z' to increase or decrease by 100, respectively");
Serial.println("'s' or 'x' to increase or decrease by 10, respectively");
Serial.println("'d' or 'c' to increase or decrease by 1, respectively");
}
void loop() {
scale.set_scale(calibration_factor); //Adjust to this calibration factor
Serial.print("Reading: ");
Serial.print(scale.get_units(), 2);
Serial.print(" grams");
Serial.print(" calibration_factor: ");
Serial.print(calibration_factor);
Serial.println();
if(Serial.available()){
char input = Serial.read();
if(input == 'a'){ calibration_factor += 100; }
else if(input == 'z'){ calibration_factor -= 100; }
else if(input == 's'){ calibration_factor += 10; }
else if(input == 'x'){ calibration_factor -= 10; }
else if(input == 'd'){ calibration_factor += 1; }
else if(input == 'c'){ calibration_factor -= 1; }
}
}
```

```
#include "HX711.h"
#include <LiquidCrystal.h>
#define CLK A0
#define DOUT A1
#define calibration_factor 2153.0
HX711 scale;
LiquidCrystal lcd(12, 11, 5, 4, 3, 2);
int resetPin = 10;
float weight;
void setup() {
lcd.begin(16, 2);
lcd.print("Initializing...");
delay(100);
pinMode(resetPin, INPUT);
scale.begin(DOUT, CLK);
scale.set_scale(calibration_factor);
scale.tare(); // reset weight on scale to 0 grams
lcd.clear();
}
void loop() {
weight = scale.get_units();
lcd.clear();
lcd.setCursor(0, 0);
lcd.print("Weight:");
lcd.setCursor(0, 1);
lcd.print(weight);
lcd.print(" grams");
if(digitalRead(resetPin) == 0){
scale.tare(); // reset weight on scale to 0 grams
}
scale.power_down(); // put the ADC in sleep mode
delay(10);
scale.power_up();
}
```

]]>In this blog post I will introduce a visually stunning effect called Marangoni bursting, explain the science behind it, and show how you can recreate this effect easily and cheaply.

Surface tension is created by the cohesive forces of neighbouring fluid particles at an interface and has a tendency to minimize surface area. A classical example of surface tension gradients is the meniscus formed at the boundary of a test tube and water. A concave meniscus is formed since the adhesive forces between the glass and water molecules at the boundary are greater than the cohesive forces between neighbouring water molecules away from the boundary. Conversely, mercury in a test tube forms a convex meniscus as the cohesive forces between the liquid are stronger than the adhesive forces at the boundary of the test tube.

Marangoni flow is a fluidic effect driven the gradient of surface tension between two fluids. As from above, surface tension describes the intensity at which of one type of liquid pulls on its neighbouring molecules. Hence if there is an interface of two liquids with different surface tensions, one liquid will pull more strongly on the other liquid. This imbalance of forces leads to Marangoni flow.

To demonstrate this visually, fill a shallow dish with water and grind pepper flakes on the surface. Using a toothpick deposit a small drop of dish soap in the center. The pepper flakes are instantly pulled outwards! This can be explained through a difference in surface tensions: water has a higher surface tension than dish soap. As such the water molecules exert a great force on its neighbours (in this case the pepper flakes). Before the dish soap is added the system is in static equilibrium: the water evenly pulls on each pepper flake in all directions. After adding dish soap there is an imbalance of forces (weaker forces pulling to the center and stronger forces pulling outwards), this imbalance drives the pepper flakes radially outwards. This is known as a Marangoni flow.

Similar results can be achieved by filling a dish with milk and several drops of food colouring then depositing a small drop of dish soap. The physics describing this system are the same as the water and pepper flakes.

Marangoni flows are fascinating in their own right, but before getting to the unique case of Marangoni bursting I must first review vapour pressure. Vapour pressure is the pressure a gaseous vapour exerts on its liquid counterpart when the system is in thermodynamic equilibrium. If two liquids (say, water and isopropyl alcohol) are left open at standard temperature and pressure, their difference in vapour pressure can explain the difference in rates of evaporation.

Isopropyl alcohol (IPA) has a much higher vapour pressure than water. Since IPA requires a higher pressure for an equilibrium between liquid and gas volumes to be maintained the liquid IPA will evaporate at a faster rate than water. A further discussion can be had, however only the fundamentals of vapour pressure are required to explain the difference in evaporation rates in the context of this experiment.

I've established the fundamental principles that can be used to explain Marangoni bursting but before explaining the physics any further, I want to provide the steps so you can reproduce this yourself.

Fill a flat, shallow dish with a cooking oil to a height of about 1cm. Separately mix water, isopropyl alcohol and food colouring, while I encourage independent experimentation a good starting point is a ratio of 3ml to 1.5ml to 2 drops of the three liquids, respectively. Use a syringe to deposit a small amount of this mixture on the bed of oil. The effect will begin instantly, sit back and enjoy the mesmerizing show.

** Tips**:

- Isopropyl alcohol evaporates rapidly, thus if you wait between mixing and depositing the mixture on the oil, the concentration of alcohol will be different to when it was initially created the mixture.
- If you are doing repeated trials, I found that to save oil you can briefly lay a paper towel on the surface to absorb the water and ethanol droplets while keeping most of the oil in the bottom of the dish for the next trail. While I didn’t test it, this could also be a unique way of making tie dye shirts.

Isopropyl alcohol has a lower surface tension than water and a higher vapour pressure, meaning it evaporates much faster. Thus, at the edge of the droplet the alcohol evaporates quickly resulting in a lower alcohol concentration at the edge relative to the centre. This lower concentration of alcohol means that the surface tension at the edge is stronger and so the interior mixture is pulled outwards via Marangoni flows. This cycle repeats with more and more of the mixture collecting at the outer rim.

As liquid collects at the outer rim, surface perturbations create instabilities and long fingers begin to protrude radially outwards. These fingers break into tiny droplets that shoot radially outwards. Together, the differences in surface tension and vapour pressure results in the stunning creation of hundreds of smaller droplets.

I encourage you to consider how the principles learnt above can be translated to explain the *“tears of wine”* commonly observed in red wine. Hint: red wine is basically just water and ethanol.

I won’t explain this effect, but with some thought and experimentation you should be able to describe the process using the principles taught above. You can check your conclusions with these great explanation videos by Applied Science and Dan Quinn

This project and blog post was inspired by the following paper:

L. Keiser, H. Bense, P. Colinet, J. Bico, and E. Reyssat. Marangoni Bursting: Evaporation-Induced Emulsification of Binary Mixtures on a Liquid Layer. PRL 118, 074504 (2017).