I’m working on a paper about an algorithm for hotstarting nonlinear program solves; one application of this might be in the realm of nonlinear model predictive control. In these types of models, we first define the physical equations for the system under consideration. They are subject to some control parameters, which are just a mathematical representation of the input we could give the system. We also define an objective–something that we would like to minimize(usually something like time or energy).

The problem we picked to test our code with is a simple one: We have several masses connected by springs, arranged left-to-right. The far left ball is fixed in position. The masses all start in a horizontal configuration. Gravity then takes over, and the masses fall into place. Since they are connected by springs, it takes a while to settle into place. We can control this system by moving the far right ball around. The goal we define is to come to rest as soon as possible(the idea being that we could move that far-right ball around to help it come to rest.)

There’s a bit more detail to the mathematical formulation, but that’s all in the paper. Here, I want to go over the visualization. So, to start, here’s the visualization I generated:

The basic flow I used to generate it: . Run the mathematical model to actually generate the optimal controls and resulting point mass positions, and output them to a file. . Run POVray on each file to generate pictures. . String the pictures together using ffmpeg.

Before moving on, I want to note that this is a really inefficient way to do this. POVray can do file reading, but I couldn’t get it working properly, and I was trying to crank something out really quickly, so I used this kinda hacky way. Hopefully someday I can update it to do a more reasonable job.

Now for the details:

Position Generation

Again, a lot of the details of the algorithm will be in the paper. The main trick I used here was to have the MATLAB code actually output POVray code. In particular, each NMPC solution generated something like:

// filename: spiral2-000.pov
#include "spiral2_common.inc"

//  1
object {Make_Spiral2(1.400397, 0.000000,0.000000,0.000000, 0.744091,-1.186355,0.000000, <0,0,0.000000> ) } 
object {Make_Spiral2(0.780335, 0.744091,-1.186355,0.000000, 1.499604,-1.381607,0.000000, <0,0,0.000000> ) } 
object {Make_Spiral2(0.750502, 1.499604,-1.381607,0.000000, 2.250000,-1.394229,0.000000, <0,0,0.000000> ) } 
object {Make_Spiral2(0.750502, 2.250000,-1.394229,0.000000, 3.000396,-1.381607,0.000000, <0,0,0.000000> ) } 
object {Make_Spiral2(0.780335, 3.000396,-1.381607,0.000000, 3.755909,-1.186355,0.000000, <0,0,0.000000> ) } 
object {Make_Spiral2(1.400397, 3.755909,-1.186355,0.000000, 4.500000,0.000000,0.000000, <0,0,0.000000> ) } 

sphere { <0.000000,0.000000,0.000000>, 0.25
        texture{ pigment{ color rgb<0,1,0>}
                 finish { phong 1}}
      }
sphere { <4.500000,0.000000,0.000000>, 0.3
        texture{ pigment{ color rgb<1,0,0>}
                 finish { phong 1}}
      }

Here, the numbers actually come from the simulation. Also note that I needed to generate some 57 different .pov files, one for each control we solved for.

POVray Spirals

So I had all the positions written out into different files, but we still need to define the Make_Spiral2 function.

I searched around and found Friedrich A. Lohmueller Spiral Tutoral. Basically, he rotates a sphere through three-space, taking the union of the many spheres, yielding a spring or spiral shape. It was pretty close to what I wanted, except that I wanted to specify start and end points for my springs and masses. So I changed the code to automatically rotate the correct amount based on the start/end coordinates given. The full file spiral2_common.inc also needs to define light sources, camera views, and any colors and textures:

// spiral2_common.inc
// POV-Ray 3.6/3.7 Scene File "spiral0.pov"
// created by Friedrich A. Lohmueller, 2003 / 2010 / Jan-2011
// Modified by traviscj, July 2013.
// Demonstrates the animation of a spiral
//--------------------------------------------------------------------------
#version 3.6; // 3.7;
global_settings{ assumed_gamma 1.0 }
#default{ finish{ ambient 0.1 diffuse 0.9 conserve_energy}}
//--------------------------------------------------------------------------
#include "colors.inc"
#include "textures.inc"
//--------------------------------------------------------------------------
// camera ------------------------------------------------------------------
#declare Camera_0 = camera {angle 25
                            location  <2.5 , -6 ,-50>
                            right     x*image_width/image_height
                            look_at   <2.5 , -6 , 0.0>}
camera{Camera_0}
// sun ---------------------------------------------------------------------
light_source{< 1500,2500,-2500> color White}
// sky ---------------------------------------------------------------------
sky_sphere { pigment { color rgb <1.0,1.0,1.0>
                     } // end of pigment
           } //end of skysphere
//--------------------------------------------------------------------------
//---------------------------- objects in scene ----------------------------
//--------------------------------------------------------------------------
#macro Make_Spiral2(Sp_Length, xx1,yy1,zz1, xx2,yy2,zz2, rotrot)
#declare Spiral  =  //--------------------------------- the spiral
union{
// #declare Sp_Length = 2.0;
 #local SPIRAL_RADIUS = .1;
 #local N_per_Rev = 500;   // Number of Elements per revolutions
 #local N_of_Rev  = 8.00;  // Total number of revolutions
 #local H_per_Ref = Sp_Length / N_of_Rev;// Height per revolution
 #local Nr = 0;                          // start loop
 #while (Nr< N_per_Rev*N_of_Rev)
   sphere{ <0,0,0>,0.025
           translate<SPIRAL_RADIUS, -Nr*H_per_Ref/N_per_Rev, 0>
           rotate<0,  Nr * 360/N_per_Rev,0>
           texture{ Chrome_Metal
                    finish { phong 1}}
         }
 #local Nr = Nr + 1;    // next Nr
 #end // --------------------------------- end of loop

  sphere { <0,0,0>, 0.2
          texture{ pigment{ color rgb<0,0,1>}
                   finish { phong 1}}
        }
} // end of union  -------------------------------- end of spiral
Spiral rotate(<0,0,90+degrees(atan2((yy2-yy1),(xx2-xx1)))>) translate(< xx1,yy1,zz1>)
#end

Running POVray

Before compiling the POVray code, obviously, POVray is required. On OSX, Homebrew handles this nicely:

brew install povray

Before running POVray, I created .ini files for each one. I think there must be a better way, but it’s a start:

// spiral2-001.ini

Width=1600
Height=1200

Antialias=On
Antialias_Threshold=0.3
Antialias_Depth=3

Input_File_Name=spiral2-001.pov

With all of those .ini files in place, the most straightforward way to run them all is just

for f in $(ls spiral2-*.ini); do povray $f; done

This gives me a folder full of PNG files.

Assembling video

This ended up being trickier than I thought it should be. I couldn’t really find a really solid video convert command, but I finally found mencoder’s image encoding syntax and ran

$ mencoder mf://\*.png -mf w=800:h=600:fps=4:type=png -ovc lavc -lavcopts vcodec=mpeg4:mbd=2:trell -oac copy -o output.avi

Note that this downsamples the resolution; I could get a bit better video by not doing that.

As usual, I uploaded everything to github.