September 8, 2021

Closed Loop Espresso Part 2: Firmware Estimation and Control

I've made a lot of progress on the software and control side of the espresso machine, and learned some good lessons about designing the real version of the machine.

This is the first post in a series about the firmware/software, more to come soon.

To "spill the beans", as people at work say, here's a demo of the current state of the machine:

Here's what's going on while making the espresso in the video:
  • After the "Start" button is pressed, water is pumped through the heater and back into the tank at a constant flow rate, and the water/group heaters start heating.
  • Once both temperatures have converged to their setpoint, flow is switched from from the tank to the group.  At first, the water is flushed through the group to the drip tray, to purge air from the group.
  • Once the air is purged, the valve to the drip tray closes and "preinfusion" starts.  The group is filled up at a constant flow rate.
  • Once the desired shot pressure is reached (indicating the group is full and puck is saturated with water), the machine switches over to pressure control, and holds a constant pressure for the remainder of the shot.
  • Once the desired output weight is reached (for now estimated by integrating a pump displacement/leak model), flow to the group is controlled to zero and excess pressure is vented to the drip tray, so the drips through the puck stop quickly.
On to the subject of this post, firmware-level estimation and control

Flow Estimation
Originally I was planning on finding a flow sensor, but so far haven't found a sensor that's both affordable and any good - the impeller type sensors for all-in-one coffee machines have pretty terrible accuracy.  Fortunately I have enough instrumentation to estimate the flow pretty well.  Here's the basic model I'm using for flow estimation:

There's an "ideal" pump with constant displacement \(D\), spinning at angular velocity \(\dot{\theta}\), with the flow characteristic:

\[F_{ideal} = \dot{\theta} D\]

There's leak path between the outlet and inlet of the pump, with a nonlinear resistance which behaves quadratically.  The pressure, \(\Delta P\) across the resistor is:

\[\Delta P = C_{1}F_{leak} + C_{2}F_{leak}^{2}\]

Solving the quadratic equation:

\[F_{leak} = \frac{-C_{1} + \sqrt{C_{1}^{2} + 4C_{2}\Delta P}}{2 C_{2}}\]

\(\Delta P\) and \(\dot{\theta}\) are both measured, and the output flow is estimated by just adding up \(F_{ideal}\) and \(F_{leak}\).

This is (of course) in SI units, so \(\dot{\theta}\) is in cubic meters, \(D\) is in cubic meters per radian, \(\Delta P\) is in Pascals, and all the flows are in cubic meters per second.

Pump Characterization

The pump displacement is listed on the datasheet, and one could probably get a decent guess at the leak coefficients by squinting at the pressure/flow curves in the datasheet, but I measured these parameters in-place.

I measured the pump displacement by just pumping water into a cup on a scale, and measuring the weight in the cup vs pump angle.  Averaging over a minute of pumping, I got a result that was nearly identical to the datasheet value (hooray!).

I measured the leak coefficients \(C_{1}\) and \(C_{2}\) by blocking off the output of the pump with the needle valve, and varying the speed of the pump.  Blocking the output sets \(F_{leak} = F_{ideal}\), so the coefficients are found by fitting a curve to \(\Delta P\) vs \(\dot{\theta}\).

Looks pretty quadratic to me:

Flow Control

Flow control falls right out of the flow equations above.  Knowing \(\Delta P\) across the pump and have a desired flow rate \(F_{des}\), the desired pump velocity \(\dot{\theta}_{des}\) is set to:
\[\dot{\theta}_{des} = \frac{F_{des} + F_{leak}}{D}\]
Then the closed-loop velocity control is handled on the motor drive.  

Pressure Control

Pressure control was a bit more involved.  I started out by measuring the pressure frequency response of the system with the puck simulator installed and blocked off for zero flow, by applying a chirp torque signal and measuring the resulting pressure

The time series data can be turned into a bode plot by taking the ratios of the FFT's of the input and output.  I've fit a 2nd order model to it, which is what I'd expect with the output blocked - the dynamics should be dominated by the inertia of the pump motor and magnetic coupling, and the spring constant of the magnetic coupling + fluid.  When the output isn't blocked off, a third pole shows up, but that didn't matter too much for designing the controller.

I designed a pressure controller by loop shaping a PI + lead controller.  I was able to push it to ~90 rad/s (14.3 hz crossover) on the hardware before things got crunchy.

Here's the closed-loop step response (compared to the simulated step response given the 2nd order fit of the dynamics and the controller above).  There's a little periodic wiggle in the pressure (I think from the teeth of the gear pump), but overall it looks pretty good.

And here's the measured vs expected closed-loop frequency response, measured the same was as the open-loop frequency response with a chirp input:

Temperature Control

As boring as temperature control is, it turned out to be a real pain.  This prototype hardware (due to being thrown together out of mostly off-the-shelf sensors/fittings/etc) has the classic problem of non-collocated actuators and sensors.  Basically, the temperature sensors are relatively far away from the heaters, so there's significant dynamics between the heaters and sensors.  That, coupled with relatively slow time constants meant testing was very time consuming.  I frequently found myself filling up the water tank with ice cubes to cool things down faster.

There were a couple mildly interesting issues I ran into though.  One had to do with the SSR's I'm using to turn on and off the heating elements.  Typical SSR's only turn on and off at mains zero-crossings.  I just stuffed some PWM on the input to the SSR, in the hopes that the probability of the input to the SSR being high at the zero crossing would be equal to the PWM duty cycle.  Depending on the PWM frequency, this could work out.  But due to some timer constraints with the clock configuration for USB and CAN, my PWM frequency ended up evenly dividing into 60 Hz.  As a result, the PWM rising had a consistent alignment with the mains zero crossing, causing the SSR to turn on much more often than I expected given the PWM duty cycle.  Rather than trying to reconfigure the timers, I implemented software PWM so I could make the PWM period much larger than 1/60th of a second, to avoid this issue.

The other interesting issue I ran into had to do with my crappy hacked-together water heater design.  Here's a cartoon of (roughly) what happens in the water heater.  Water flows in at some inlet temperature, and out (hopefully) at the desired temperature.  Assuming the heater itself is a constant temperature, which is probably not an unreasonable assumption if it's highly thermally conductive, then the average water outlet temperature exponentially approaches the heater temperature throughout the heater.

Average is the key word.  As you can see from the cartoon plot, depending on the length of the heater and the convection coefficient between the heater wall and the water, the heater temperature could potentially be much higher than the outlet temperature.  And in reality, at a given cross section of the water flowing through the heater, the entire cross section will not be at the average temperature.  The water towards the walls of the heater will be hotter than the average, and the water in the middle of the cross section will be cooler than average.

For making espresso, the desired output temperature is typically somewhere in the low-90's C, which is awfully close to water's boiling temperature under ambient pressure.  The larger the difference between the heater temperature and the average outlet temperature, the larger the temperature gradient across the cross section of the water flow - so with a bad heater design (ahem) where the heater needs to be significantly hotter than the outlet, the water at the perimeter of the flow cross section can actually start boiling, even though the average water temperature is below boiling.  This results in the machine spluttering as it cycles water back to the tank while preheating, and introduces (highly compressible) steam into the group when the shot is happening, which affects the pressure control.

This issues seems solvable, and I think with a fresh heater and group design I can also significantly decrease the thermal mass and improve the temperature control bandwidth.  

This is getting long, so I'll pick up next time with a software post.

April 15, 2021

New motor drive firmware

 Over a year ago I started porting my motor control firmware away from mbed.  I've been working on it in spurts every few months and and ended up re-doing a lot of things from scratch, so it took much longer than I originally planned.  Find it here:

The core motor control math hasn't changed (foc.c is nearly identical to the old version).  Lots of features have been improved though:

  • Encoder calibration/linearization is now current-mode, not voltage mode, and the calibration current can be configured from the serial terminal interface
  • Lots more parameters are configurable from the serial interface
  • Hardware setup is done with the Cube auto-generated setup code.  Linking between the hardware and motor control code is all done in hw_config.h so it should be pretty easy to port to different micros (someday I'll switch to using G4's instead of F446's....)
It's certainly not completely polished, but I think everything's working - I've been using this firmware for the espresso machine pump.  I don't have any Mini Cheetah hardware other than the motor drives any more, so I haven't been able to confirm backwards-compatibility.  Let me know if you try it and have problems.

February 25, 2021

Closed Loop Espresso Intro

Even before getting a lever espresso machine, I've been thinking about how I might design one if I were to do it from scratch.

This is the layout I've converged to:

The basic summary is:
- Electric motor driven pump pumps water through a heater.
- If the water at the output of the heater isn't up to temperature, the 3-way valve after the heater will cycle the water back through the tank.
- Once the water at the output of the heater is up to temperature, the valve switches and water is piped over to the group head and coffee.
- An extra "drain" valve between the group head and the drip tray relieves pressure after the shot is done.  Strictly this isn't necessary, but it serves a few purposes I'll get into later.

Critically with this setup, the whole tank of water doesn't have to get hot as it would in a boiler-based espresso machine.  I'm hoping to get the delay between powering the machine on and pulling the first shot down to a minute or less.  With a standard 15A mains outlet, it should be possible to real-time heat over 5 mL/s of water from room-temperature to espresso temperature, so main delay is in getting the thermal mass of the heater up to temperature.

The pump:

As far as I can tell, most "real" espresso machines use either vibratory pumps or rotary vane pumps.  From a pure performance perspective these don't seem like great solutions to me.  Vibratory pumps are noisy and hard to control (not saying you can't, Decent does pressure control with vibratory pump), and their main appeal is they're cheap and you can plug them straight into mains.  Rotary vane pumps are positive displacement (which is good for pressure control, roughly speaking), but tend to have large displacements - meaning for a given pressure, a lot of torque is needed at the input, since (ignoring losses) \(energy = torque*angle = pressure*displacement\).  Lots of torque = large motor or gear reduction.  However, the pumping power requirements for making espresso are tiny - on the order of a few watts of \(pressure*flow rate\) (9 Bar * 4 mL/s is 3.6 watts, and I think that's towards the upper end of espresso flow rates).  So a pump with a small displacement and a tiny motor (effectively getting "gear reduction" through the pump) makes more sense to me.  

Gear pumps were my first though.  What I wanted was "stainless steel body gear pump with PEEK gears", and it turns out that exactly this object exists and can be purchased on ebay for a few tens of dollars.  Even better, they all are magnetically coupled, so there are no shaft seals to leak.  A couple weeks later I had a Fluid-o-Tech MG304 pump.

Here's a shot of the inside.  The gears and bushings are CF-PEEK and the rest is stainless:

Here's a relevant plot from the datasheet.  The blue curve labeled "A" is for a 1 cP viscosity fluid, e.g. water.    Should be able to get over 200 psi at greater-than-espresso flow rates.  Since the pump displacement is so small, torques needed are in the ~.1 N-m or less range, which means I can use a tiny motor.

The main downside of gear pumps is leakage - which you can see in the pressure vs flow plot.  Gear pumps "seal" by having very close tolerances between the tooth tips and faces of the gears and the pump cavity, so some fluid leaks through these gaps.  For this pump it's very linear with pressure - at a constant input speed (which would be a constant flow if there were no leakage), output flow drops off linearly with pressure drop across the pump. 

The pump came with a low-end brushed maxon driving it.  I swapped that for a Moog BN23 Silencer motor I snagged from a Media Lab trash pile ~7 years ago.  One of its larger brethren was used to power the tiny lathe, and another by Dane for Robot Art.  A 3d-printed spacer adapts it to the pump.  Eventually this'll need to be replaced with a more temperature resistant part, since the water flowing through the pump will get fairly warm.

I modified a spare pendulum motor controller so I could plug in the pressure sensor.  Eventually I'll  have a separate micro for doing all the high-level control and use one of my usual motor drives with the built in encoder, but for now this was the easiest thing to hack together.  

This was my pressure control test bed -  the pump and motor, an automotive-style thread-in pressure sensor, a dial pressure gauge, and a needle valve as a variable load;

Testing closed-loop pressure control.  The pressure is plotted/setpoint set from the computer, and you can see the response to the needle valve opening and closing.  With this setup I can get around 10 Hz of closed-loop pressure bandwidth.

Some kind of interesting notes about the pressure control, looping back to the pressure relief valve.  One of the factors limiting the pressure control bandwidth is the compliance in the pump system.  A lot of this is from the magnetic coupling, which has a sinusoidal force-vs-displacement profile, with zero stiffness around zero displacement, and relatively low peak stiffness.  The other source of compliance is the fluid circuit - the bulk modulus of the fluid, the stiffness of the tubing and other stuff between the pump and the output, and much more importantly, the volume of air trapped in the system.  Without the relief valve, a big air bubble could get trapped in the top of the group head, above the puck of coffee, which would add a huge compliance.  To avoid this, the relief valve can be kept open until the group head is full of water, and then closed.

On to water heating:

The water heater is made from three 1/4" x 7" die heaters sandwiched between two thin aluminum plates, with a stainless tube snaking between the heaters.  My main goal was minimizing the thermal mass of the heater.  This design is kind of terrible, but it's functional enough to get me going.  The bend radius is super tight, and I did a few test bends that went well, but most of the bends in the final heater are awful.   Still, they didn't crack and still pass water, so good enough.

Two thick sheets of FR4 sandwich the whole assembly for insulation.  The machinist clamps are for extra pressure in the center.  Yeah, this design was bad.

For water temperature control, I have a platinum RTD in the flow of water, and a thermistor measuring the temperature of the heater.  I circulate water through the heating loop and back into the tank (in the picture above the metal pot...) until the the water temperature sensor has converged to the desired temperature, and then a 3-way valve is switch to pipe the hot water to the group head instead of back to the tank.  It takes around a minute to converge from room temperature to ~90C, and that's with a pretty terribly designed heater and a water temperature sensor with way too much thermal mass (I used a generic NPT threaded sensor so I could throw this together quickly).  So I think with some actual design I'll be able to get it even faster.

Similar to the heater, the group head was designed to be low thermal mass.  This is the opposite of how most espresso machines work, where the groups often have as much thermal mass as possible for "thermal stability"  But they're usually heated by water or steam from the boiler.  With the power of feedback control, low thermal mass, a group heater, and a group temperature sensor seems like the way to go.

I designed the group to work with a La Pavoni  portafilter, shower screen, and group seal, since I already had those lying around.  For the final version of all this (assuming there is a final version), I'll probably switch over to the standard 58mm portafilter size.

Here's the CAD and a cross section.  The typical wall thickness is 1.5-2mm.  I could definitely have gone harder with the light-weighting, but I wasn't in the mood for any analysis (this is pressurized to ~150 psi), so I just eyeballed it for now.  Next round of parts I'll probably be able to make it ~50% lighter. 

A screenshot of the CAM for making the top half of the group.  Thanks to the power of the 5-axis mill, it was only 2 setups, despite the holes pointed in every direction.

Video of the machining:

The finished parts assembled with seal and shower screen;

With the group head done, all it took was an afternoon of plumbing, sketchy wiring, and a tiny bit of firmware to read a few new temperature sensors and toggle valves.

Everything's just bolted to a strip of aluminum right now, and I stole the drip tray from the La Pavoni.

Here's the first test of the machine.  In my excitement I let it run for a little too long and the shot ended up over-extracted, but it turned out way better than I expected for a first try.  Pardon the wrench used as a handle for the group head - McMaster order with with a threaded handle was delayed.

For this test the group head wasn't heated, but I do have a 100W cartridge heater and thermistor for it, so the water flowing through it won't be cooled down.

While right now the machine is a huge heap on a benchtop, I think with some actual design it could be packaged to be very compact.  

December 28, 2020

New Photo and Video Feed

 You'll now notice that I added a Feed tab to the top of the page.

When I started the blog, I would wait until I completely finished a project, then do one gigantic writeup at the end.  This habit was left over from when I used Instructables to document things.  As my projects got more complicated, that blogging style stopped working.   So I started the In-Progress section for shorter technical updates. When I (on occasion) finish something, I now throw a summary on the main page, and links to the relevant posts from the in-progress section.  I think this has been working fine, but I still tend to hold onto stuff until I've made enough progress for a proper technical blog post.  

Despite the relative infrequency at which I write detailed posts, I think most weeks I do something side-project related that someone on the internet would find interesting.  To avoid cluttering the proper blog, I made a new Feed section where I plan on just throwing up images or videos with one or two sentence descriptions.  "Hey, that sounds like Instagram", you might say, and you'd be right.  But I've avoided Facebook and Instagram this long, so I don't intend to stop now.  Maybe I'll get bored of this and it will fall by the wayside, but for now, enjoy.

December 9, 2020

La Pavoni Europiccola Lever Espresso Machine Restoration

I've had a La Pavoni Professional lever espresso machine for a couple years now, which I got on ebay and fixed up.  Getting it working turned out to be completely uneventful, so I never wrote up anything about it.  I recently restored a Europiccola for my sister, and this time around the restoration was much more involved.  It was in pretty rough shape when I got it - all the seals were shot, the sightglass was smashed, and the copper/brass plating was  peeling horribly.  But the heating element and electronics were in good shape, and core mechanical pieces worked.

I forgot to fully document the state I received it in, but here's what it looked like minus the grouphead:

All the seals needed replacing.  The trickiest of these was the seal in the top of the grouphead.  This seal is held in by a brass washer and a retaining retaining ring, which was made out of non-stainless steel.  It sits in a chamber full of steam, so I don't know what they were thinking there.  So the retaining ring was extremely rusted, so much so that the holes used for removing it with snap ring pliers were gone:

I was able to slowly chisel it loose and pry it out with a dental pick.  This is what was left of the original retaining ring (left), seal, and washer:

Here it is after installing a new seal, cleaning up the brass washer, and using a new stainless retaining ring.  Sorry for the blurry pictures:

This machine originally had a copper finish on the boiler and gold/brass finish on the grouphead and base, but both were discolored and peeling.  The finish on this  unit was strange -  on the outside, there was a clear varnish layer (the source of the peeling).  Under that, the boiler has an electroplated layer of actual copper, and the base has an electroplated brass layer.  Neither of those layers were in very good shape.  Where the copper/brass were wearing through, there were signs of a chrome coat beneath that. 

I took the boiler as far apart as I could before stripping the coating.  The sightglass attachment points are bolted through the boiler, with their nuts on the inside of the boiler, at a funny angle.  I'm sure there's a specialized wrench for removing these, but I didn't have that.  I used a 16mm wrench with a 12-point ratchet.  Because the opening at the base of the boiler is pretty small, the wrench couldn't be turned enough to click past one pawl of the ratchet.  To make progress, I had to turn the wrench a few degrees, remove the wrench, manually click the ratchet over two clicks, take up the backlash in the ratchet, put it back in, and turn the nut a few degrees.  And repeat several dozen times.  I couldn't get the wrench around the nuts holding on the pressure relief valve or steam wand, so I just gave up and left them on.

The original plated bolts holding that hold the grouphead on were a little rusty, so I got some new stainless steel bolts, turned the text off the caps, and gave them a polish:

A combination many rounds of paint stripper and scrubbing with a gentle polishing compound stripped off the varnish and electroplated layers, and revealed an almost pristine chrome coat on the boiler.  The base wasn't in quite as good shape, but a few cycles of scrubbing with scotch-brite and sanding with progressively finer sandpapers improved it.

Part way through stripping, the boiler had this really interesting crackle-texture:

The chrome on the base slowly being revealed:

The drip tray had a few deep rust pits in it.  I stripped the rust out of the pits, and sprayed a layer of clear high-temperature spray paint to over them hopefully keep them from eating all the way through the base.  Fun fact I learned in this process, some "high strength" toilet bowl cleaners are 5-10% hydrochloric acid, so it works really well for stripping rust or black oxide from steel.  

Here's what the pits looked like after a couple rounds cleaner:

I was amazed at the state of the chrome coat underneath all the gunk.  Going in, I was sure I would need to strip the parts all the way down and repaint them to make it look decent.  Lots of scrubbing later, here's what the boiler looked like:

Kind of interesting side note - I could have taken off all the copper with just paint stripper.  I noticed that when I wiped off the paint stripper residue, it was blue from oxidizing the copper.  Scrubbing was much faster though.  I used some glass stovetop cleaner, which is very gently abrasive and won't scratch the chrome.  Steel wool probably would have been much faster, but I didn't have any on hand.  Steel is softer than chrome, so I don't think it would have scratched the chrome.  

The base wasn't in quite as good shape as the boiler, but was able to get most of the blemishes out with scotch-brite and sandpaper.  In the end it doesn't look as nice as the glossy boiler, but the smooth brushed look is still much better than it started out.  Here it is going back together:

I stuck the logo back on with some 3M VHB foam tape:

The fully assembled machine:

Once it was back together I ran two rounds of Dezcal descaling solution through it, did two rounds of letting it sit full water with a couple spoonfuls of baking soda dissolved in,  and filled/boiled/flushed the water probably a dozen times to make sure dirt and all the chemicals I used were out of the system. 

Then it was ready to go.  Here's pulling a shot of Red Bird Espresso to test it out:

November 17, 2020

Transmission Ratio Trajectory Optimization

 Getting back to this mechanism.  In the post about the dynamics with varying transmission ratio, I came up with the design for the transmission profile by a combination of hand-calculations, intuition, and guesswork.  But we can do better.  

This problem fits nicely into a trajectory optimization.  For an intro to trajectory optimization, I'd highly recommend this MIT Underactuated Robotics chapter

I not-so-subtly hinted that the point of this mechanism was for jumping or launching things with an electric motor.  The simple English explanation of what I'm trying to do with a trajectory optimization here is answer this question:  What transmission ratio vs time maximizes jump height without breaking the mechanism? 

The optimization-lingo version of that would be something like:

Minimize \(-\dot{x}(t_{final})\), subject to:

  • Dynamics are not violated (i.e. \(x(t) = \int_{t_{0}}^{t_{final}} \dot{x}\,dt + x(t_{0})\))
  • Initial conditions (everything starts at rest, i.e. \(x(t_{0})=0\), \(\dot{x}(t_{0})=0\), etc)
  • \(\ddot{x} < a_{max}\), \(\omega < \omega_{max}\) (acceleration limit to limit forces, maximum motor speed)
  • And a few other things I'll get to later

I set up the problem using CasAdi, which is the same tool Jared and I used for the Mini Cheetah backflip optimization.  CasAdi makes it super easy to set up nonlinear optimizations like without having to understand what's going on in the back-end too deeply.  

Walking through my code (on GitHub here):

Pick some constants

The only "choice" here is the number of trajectory intervals.  Roughly speaking, too few will result in integration error, and too many will take longer to converge.  Everything else is just physical parameters:


N = 200              # number of trajectory intervals

m = .75              # mass            (kg)
j_rotor = 6e-6 	     # rotor inertia   (kg*m^2)
l_leg = .35 	     # leg length      (m)
tau = 1.6 	     # motor torque    (N-m)
w_max = 1200         # max motor speed (rad/s)
g = 9.8 	     # gravity         (m/s^2)

Decision variables

I took a "Direct Transcription" approach to the optimization - at each point along the trajectory, both the state (positions, velocities) and control inputs (transmission ratio and its derivative) are decision variables being optimized for.  As someone with not very much experience with this stuff, I found this easier to understand and implement than other methods (see here for several formulations).  I end up with a ton of redundant decision variables here (really there's only one independent one, which is the transmission ratio), but having position and velocity and acceleration and motor angular velocity and transmission ratio derivative as decision variables makes it very easy to put constraints on those terms.

Another note here, I have the final takeoff time as a decision variable.  This seems to usually not be a good idea (I won't pretend I can do an explanation justice), but this problem is simple enough that it works out.  Doing the Mini Cheetah backflips, for example, we did not do this, and fixed the takeoff times ahead of time (and did a little manual tuning to get the optimization to produce a nice result).  But that optimization had way more states and control inputs than this one.

Also, CasAdi is great.  foo = opti.variable(rows, cols), and you have a matrix of decision variables.

opti = Opti() 		# Optimization problem

###  decision variables   ###
X = opti.variable(2,N+1) # position and velocity trajectory
pos   = X[0,:]
speed = X[1,:]

accel = opti.variable(1, N+1) # separating this out behaved better than including in X

U = opti.variable(2,N+1)   # transmission ratio and its derivative
k = U[0,:]
kd = U[1,:]

W = opti.variable(2, N+1)  # rotor angle and angular velocity
theta = W[0,:]
thetad = W[1,:]

T = opti.variable()      # final time

Pick a cost function

Often it's hard to describe what you "want" the optimization to do in an equation, have multiple competing objectives, and so on, but this problem is charmingly simple.  Maximizing jump height is the same as maximizing takeoff velocity, so we just minimize the negative of that.

#### objective  ###
opti.minimize(-(speed[-1])) # maximize speed at takeoff

Set up the dynamics constraints

This is the real meat of the optimization.  It implements the varying transmission ratio dynamics (with gravity this time) at every timestep along the trajectory.

The equations of motion,

\ddot{x} = \frac{k\tau + j\dot{k}\dot{\theta}-k^{2}mg}{j + k^{2}m}

are implemented by the function \(f(x, u)\).  \(x\) is really the vector \([x, \dot{x}]\) at the timestep being evaluated, and u, the control input, is the transmission ratio and its derivative, \([k, \dot{k}]\).  It returns \([\dot{x}, \ddot{x}]\).

To ensure the dynamics are respected, there's trapezoidal integration constraints at each timestep, for all the states.  Basically, \(x(t+\Delta t) = x(t) + \Delta t\frac{(\dot{x}(t) + \dot{x}(t+\Delta t)}{2}\), and same deal for \(\dot{x}\), \(k\), \(\dot{k}\), and \(\omega\).

Finally there's the transmission ratio constraint, \(\dot{x} = k\dot{\theta}\)

####  dynamic constraints   ###
f = lambda x,u: vertcat(x[1], (u[0]*tau + j_rotor*x[1]*u[1]/u[0] - m*g*(u[0]**2))/(j_rotor + m_body*(u[0]**2))) # dx/dt = f(x,u)

dt = T/N # length of a control intervals
for i in range(N): # loop over control intervals
   k1 = f(X[:,i],  U[:,i])
   k2 = f(X[:,i+1], U[:,i+1])
   x_next = X[:,i] + dt*(k1+k2)/2            # trapezoidal integration
   opti.subject_to(X[:,i+1]==x_next)         # dynamics integration constraint
   opti.subject_to(accel[i] == k1[1])			# acceleration

   k_next = k[i] + dt*(kd[i] + kd[i+1])/2	   # transmission ratio integration constraint

   theta_next = theta[i] + dt*(thetad[i] + thetad[i+1])/2 
   opti.subject_to(theta[i+1] == theta_next)	   # angle/angular velocity integration constraint
   opti.subject_to(thetad[i] == speed[i]/k[i])	# linear and angular velocity transmission ratio constraint

Bounds and boundary conditions:

The bounds are hard limits on the range of values the decision variables are allowed to take.  These are mostly based on physical constraints of building a variable transmission.  There are bounds on the derivative of the transmission ratio, so it doesn't instantaneously jump.  There are bounds on the transmission ratio itself, because it can't be infinite or too small.  The motor has a maximum angular velocity.  There's an acceleration limit to limit the forces on the transmission.

In addition to the bounds, there are some constraints on the initial and final state.  Everything has to start at rest at 0 position and velocity.  At the end, the extension of the "leg" should be the length of the leg.  The final motor speed should be small, so it doesn't crash into a hard-stop with too much energy.

### bounds  ###
opti.subject_to(opti.bounded(0,kd,50))             # transmission ratio derivative bounds (no infinite slopes)
opti.subject_to(opti.bounded(0.0008,k,1))           # transmission ratio bounds (meters per radian)
opti.subject_to(opti.bounded(0, thetad, w_max))	   # maximum motor angular velocity
opti.subject_to(opti.bounded(.01, T, .5))          # reasonable takeoff time limits
opti.subject_to(opti.bounded(0, accel, 1000))      # acceleration limits

####  boundary conditions  ###
opti.subject_to(pos[0]==0)       # 0 initial position
opti.subject_to(speed[0]==0)     # 0 initial velocity
opti.subject_to(pos[-1]==l_leg)  # final position = leg length at takeoff
opti.subject_to(theta[0] == 0)	# initial rotor angle of zero
opti.subject_to(thetad[-1]<200)	# final rotor speed < 200 rad/s
opti.subject_to(thetad[-1] == speed[-1]/k[-1])	# linear and angular velocity transmission ratio constraint at the end
#### misc. constraints   ###
opti.subject_to(T>=0) # Time must be positive

Initial guess

The initial guess is the values for the decision variables that the optimization algorithm starts from.  For weird nonlinear trajectory optimizations, good initial guesses can make an enormous difference in the convergence speed and quality of the solutions.  A good example from my old lab is Gerardo's regularized predictive controller, which, as a flavor non-linear model predictive controller (MPC), runs a trajectory optimization at ~100 hz, so it can be used for (incredibly impressive) feedback control.  One of the many things that lets that optimization run at 100 hz is a good initial guess (a.k.a. "warm starting").

But this problem's simple, so I did a few runs starting from guesses of zero for everything, and picked what looked like the "average" values of the solutions I got as the new initial guesses.

####  initial values   ###
opti.set_initial(k, .02)
opti.set_initial(kd, 0)
opti.set_initial(T, .2)


Set up and solve the problem with the solver of your choice.  I used IPOPT because that's what we used for Cheetah stuff.  I played around with a few others and didn't get any better results.

### solve  ###
opti.solver("ipopt") # set numerical backend
sol = opti.solve()   # actual solve

And that's it.

Here an example solution:

Here you can see how the kinetic energy of the mass, of the rotor inertia, and the potential energy of the mass trade off over the trajectory:

Compared to the hand-designed version of the transmission profile, the optimized one can exactly max out my acceleration constraint, exactly max out the motor speed constraint, and rapidly decelerate the motor to my threshold at the end (ignore the acceleration at the last timestep, it's not meaningful):

Below is an example of an optimized transmission ratio vs rotor angle on the left, vs my hand-designed one on the right.  The mass/torque/inertia/stroke used for these are different, so the absolute numbers shouldn't be compared, but qualitatively, I was definitely on the right track with my hand-design.   Which is always satisfying to see:

The numbers I'm seeing out of the optimization are nearly 10 meters (if I can hit the weight target in there), so things are going to get exciting.

November 15, 2020

Music Server

Until Google axe-murdered it recently, I used Google Play Music to store my music collection and stream it to my various devices.  I didn't pay for a subscription, just uploaded all my music files to Play, and treated it as free storage, organization, and streaming.  Rather than switch over to YouTube Music (I tried, it was utter garbage) or some other streaming service, I took this opportunity to set up a personal music server.

My music buying and listening habits are a little unusual:

  • I want a digital copy of my music backed up on storage I own
  • I purchase and download full albums (not single songs)
  • I almost always listen to full albums (no shuffling, auto-generated playlists, etc).   
So I'm not really into the idea of subscription-based streaming services.  Compounding my dislike of streaming subscriptions, I want my money to actually go to the musicians who make the music I listen to, and streaming just doesn't do that very well.  

For the last few years, I've been tracking my album purchases as well as how many total song plays I have.  Here's what my cumulative music spending looks like since April 2018:

Since then, I've spent an average of $13.64 per month on buying albums, which is slightly more than a Spotify membership.  Looking at the sources I've purchased the albums from (directly from artist websites, Bandcamp, Amazon, etc), on average 76% of my spending goes towards the rights holders - $332.98 since I started keeping track.  If I had used Spotify, on the other hand, the total payout to rights holders would have been on the order of $65.60, given how many plays I've logged in that time (~20,500 plays, assuming .32 cents/stream).

So I spend 36% more a month on music than a Spotify subscription would cost, but 5 times as much of my money makes it to the music's rights holders.  Some other service than Spotify might be better (according to the source I linked to, Amazon pays out ~3.7x what Spotify does per-stream), but either way, given my listening habits, buying albums seems like the way to go from an artist-profit perspective.  

The Setup:

I'm now hosting my music on a Raspberry Pi running Airsonic.  For storage I'm using a pair of USB flash drives in mirrored RAID.  

I'm not going to write up a full how-to, but here's a list of resources I used and notes from setting it up.  It was surprisingly straight-forward, even though I'm no Linux wizard and have never self-hosted a website before.
  • Setting up the flash drive RAID array
  • Installing Airsonic on a Raspberry Pi.  The only thing I did differently from this guide was use open-jdk rather than oracle-jdk.  Skip past the "Set up a reverse proxy" portion for now.
  • I got a domain name on Google Domains, installed ddclient on the Pi, and followed these instructions to set up dynamic DNS so the domain points towards my router's dynamic IP address.
  • I struggled for a while with getting a TLS certificate using certbot.  It turns out my ISP blocks port 80, so I followed these instructions to get a certificate using port 443.
  • Airsonic Apache reverse proxy setup.  Once I could reach my domain name from outside of my local network, I followed these instructions.  This just worked, so my Airsonic login showed up at <my_domain_name>.com/airsonic.
  • Airsonic works best with an Artist/Album/SongName folder structure, so I used MP3TAG to re-organize my files, following the instruction here.  It was fast and worked great.
And that's pretty much it.  I've had no problems with the browser-based playback, and the Android app Subsonic seems like it works, but I haven't tested it much.