Anne Marsan: Software for Engineers

Advances may come from a single individual -- often young, often working alone, often with limited resources, often not at what's considered a premier institution -- who perseveres in the face of the conventional wisdom. It's the idea that counts. (Meyers, M. A.; American Journal of Roentgenology, March 2008, pp. 561-564.)

Awhile back I wrote part 1 of a series on design optimization. At long last here is Part 2.

As a quick review, design optimization is a mathematical process that we can use to (i) describe a system in terms of mathematical variables, (ii) change the values of the variables in order to explore different system configurations, and (iii) find the values of the variables that give us the "best" system. Design optimization can be applied to many different kinds of systems. If you can model your system with mathematical equations, then most likely you can apply design optimization to it. For instance, mathematical equations can be used to model the flow of a fluid such as water in a series of interconnected pipes. Design optimization can be used to find the best size for each pipe in order to minimize the cost of the pipes, but still move the the fluid through the system without backups or ruptures.

In part 1 of this series, I introduced the key concepts in modeling a system: our system is defined in terms of variables, the system must satisfy certain constraints, and we are trying to minimize an objective function. Usually design optimization problems are written like this:

Note that we group the design variables into a vector x. One thing I didn't mention in Part 1 about the design variables is that they are generally assumed to be real-valued, which is indicated in the last line of figure above.

Once we've described our optimization problem using mathematical equations, the question becomes, "How do we find the value of x that results in the best design?" The basic idea is that we start with some arbitrary value, called x_i, and then we use a systematic way of tweaking that value to get x_i+1. Ideally, our objective function evaluated at x_i+1 will be less than at x_i, or f(x_i+1) < f(x_i). Following the design algorithm, we would have a series of values x₁, x₂, ... , x_n, and corresponding series f(x₁) > f(x₂) > ... > f(x_n), with x_n being the optimal point. Schematically, it would look like this:

At this point, two questions must be answered: (1) What goes on in the "Do some math" box? And (2) how do we know when to stop iterating?

I'll answer the second question first. There are several different stopping criteria that are used in design optimization practice. One of the most common is to stop if the the objective function doesn't change very much between x_i and x_i+1, i.e. if:

This works if as we approach the optimal point we take smaller and smaller steps. Another approach is to stop if we reach some predetermined maximum number of iterations, i.e. if:

n >= max_iterations.

There are other criteria, especially depending on the optimization algorithm that is used, but the ones shown above are fairly common.

But back to question 1, what goes on inside the black box? This is where most of the work in design optimization takes place: designing and tweaking the optimization algorithm. The type of algorithm that is used depends on the type of design variables - whether they are continuous or discrete; the type of constraints - whether they are linear or non-linear; and the objective function - can derivatives be computed or not.

To give you an idea of the variety of optimization algorithms that are available to choose from, take a look at Wikipedia's summary on the topic. It lists 89 different algorithms. Some are fairly exotic and used in highly specialized situations. Others are more general purpose and are the workhorses of the design optimization community. In my next installment, I'll talk about some of the commonly used algorithms, going over the basic idea behind them and in what scenarios they can be applied.

I use the charting library TeeChart from Steema for one of my clients. This library is large and rambling, uses non-conventional naming schemes, and lacks good documentation. On the other hand, it does allow the user to do some very nice things with plots, if she can figure out how to do them. It wouldn't be my first choice for a plotting library, but the client has it and has asked me to use it, so use it I will.

In TeeChart, one adds to data to a chart in the form of a series. There are all sorts of useful, predefined series to choose from, but the most basic is a Line series. I add (x,y) points to the Line series, and when plotted, those points are connected with a straight line. By default, TeeChart sorts the points on the x-value, which for time series data makes sense. But for a lot of my work, I'm plotting generalized 2D paths, so it is imperative that the points maintain the ordering that I assign to them. To turn off sorting for a TeeChart Line series, use the following:

Steema.TeeChart.Styles.Line dataSeries = new Steema.TeeChart.Styles.Lines();

dataSeries.XValues.Order = Steema.TeeChart.Styles.ValueListOrder.None;
dataSeries.YValues.Order = Steema.TeeChart.Styles.ValueListOrder.None;

If you want your points to automatically be sorted on the y-value, then use the HorizLine series.

That problem I was able to solve after searching on Google and the Steema forums. The next two problems have been trickier.

First, I want to be able to reposition an axis Title. ("Title" is the Steema terminology; to me it should be called an axis label, but they use "Label" to mean the labels associated with axis tick marks.) By default the axis title is centered on the axis. I'd like to be able to move it along the axis because if I've positioned my axis so that it is not at the edge of the plot (e.g. if I have the axis at the center of the plot), then the data that I'm plotting might sit right on top of the axis label. According to Steema, Titles cannot be repositioned, so I'm supposed to use Annotations and work out their position based on the position of the axis and the size of its tick marks. Seems like a big pain, if you ask me!

My second problem is the color that is automatically assigned to a series. TeeChart has 11 built in color palettes, and to switch between them you do the following:

Steema.TeeChart.TChart chart;
int paletteIndex;

chart.Aspect.ColorPaletteIndex = paletteIndex;

I've also seen it done this way:

Steema.TeeChart.Themes.ColorPalettes.ApplyPalette(chart, paletteIndex);

Then when series are added to the chart, each one is assigned the next color in the palette (there seems to be about 15-20 colors in each palette). But when you get to the end of the palette, every new series gets the last color in the palette. The indexer doesn't seem to wrap around back to the start of palette. To work around this, I wrote my own code to rotate through a list of colors, and I manually set the series color each time I add a new series, using code like this:

dataSeries.Color = GetColor();

static int colorIndex = 0;

private static System.Drawing.Color GetColor()
{
    ++colorIndex;
    if (colorIndex > 3)
        colorIndex = 0;

    switch (colorIndex)
    {
        case 0:
            return System.Drawing.Color.Red;
        case 1:
            return System.Drawing.Color.Blue;
        case 2:
            return System.Drawing.Color.Green;
    }

    return System.Drawing.Color.Black;
}

This solution is fine, but seems like it should be unnecessary!

Currently I'm working on three different projects, and the one common component of all three is that they involve design optimization. Design optimization can be used in any discipline to find the "best" design for a particular system (the definition of "best" will vary from discipline to discipline, from system to system, and from one situation to another). Methodical, mathematic approaches to design optimization have been developed over the past 50 years, and their practical implementation, especially for large problems, has been made possible through the use of computers. A basic understanding of optimization is a good thing for software developers to have, especially those who are writing scientific and technical tools. This post is part one of a short primer on design optimization. I perused the book Principles of Optimal Design: Modeling and Computation by Panos Papalambros and Douglass Wilde before writing this, so their influence cannot be escaped.

Any system that can be modeled mathematically can be optimized. This means that disciplines as varied as engineering, physics, chemistry, biology, statistics, economics, finance, and others can make use of design optimization. Within engineering, systems such as automotive engines, building structures, HVAC systems, integrated circuits, etc. can be optimized. Two very simple systems that could be optimized are shown below. The first is a collection of 2D data points that are to be described by quadratic polynomial (a parabola). The second is a cantilevered beam that has a load applied at its tip. Different aspects of the beam can be described with mathematical equations, such as the geometry of the beam, the loading along the length of the beam, and the stresses in the beam at an point within it.

The next requirement for a design optimization problem is to identify a set of variables that describe design alternatives. These are parameters that the designer has control over, things that can be changed during the design process. For our data fitting problem above, the design variables would be the coefficients in the polynomial; the engineer can tweak their values to modify how the polynomial fits the data point. For the cantilevered beam problem, the designer has control over the length, height, and width of the beam (additionally, the beam material could be a design variable, but for this example we'll say that the material choice has been predetermined and the designer cannot change it).

Now we need to choose an objective function, expressed in terms of the design variables, that should be minimized. The objective function is frequently called a cost function. Note that if you ask any manager of a design optimization project what it is they are trying to minimize, they will say cost! It's up to the designers (and managers) to figure out what contributes to system cost and express it mathematically. In our beam example, clearly the larger the beam is, the more it will cost. But the beam has to perform a function, such as support a floor in a building, and therefore we can't minimize the size of the beam. Instead, a more typical objective would be to minimize beam deflection at its tip (if this beam is indeed supporting a floor in a building, it will make the building occupants feel a lot better if the floor doesn't sag when they walk over it). For our data fitting problem, the objective is to find the polynomial that best fits the data. A typical way of measuring the fit is to sum the squares of the distance of the given data points to the polynomial.

Designers next need to identify a set of constraints, expressed in terms of the design variables, that must be satisfied in order for their system to be valid. These can be practical limits on design variables, or more complicated relationships that describe the function of the system. For instance, in our data fitting problem, perhaps we want to ensure that the peak of our parabola is at x greater than some value d. Typical constraints for a beam design problem might be that the material properties are fixed because a certain material has been selected, or the stresses need to be below the limits of the material.

Now we have a complete specification for our design optimization problem: variables, objective function, and constraints. The formal way of expressing an optimization problem is:

minimize f(x)
s.t. (subject to)
g(x) < 0
h(x) = 0

The final step is to solve the problem, which is of course easier said than done! I will discuss this in the next installment.

But before I close, I want to point out that modeling a system using mathematical equations can itself be a challenging problem. In engineering, we spend much of our bachelors degree education learning how to do this, and it's something that doesn't always come easily, even to seasoned practitioners.

Since I'm new to this blogging thing, it has taken me awhile to figure out what tools to use to make writing and publishing posts easy to do. After reading what others have to say on the topic, and talking to my co-workers at SRT, I have hit upon Windows Live Writer with the Code Snippet plugin. But it hasn't been straightforward getting these tools to cooperate with Community Server 2007.1, the server SRT uses for blogs. Here are two tricks to make things work right.

Windows Live Writer makes writing a post and including things such as images and videos and bits of code very easy. One of the big advantages is that it can go to your blog site and grab your template, so that as you compose your post you can see just what it will look like when it is published. I could have used that feature when I wrote my last post on ZedGraph and made the images too wide for the main text column. To get the template, WLW creates a temporary post, publishes it to the server, then goes to that post and downloads the CSS. Seems straightforward...except that it failed the first and second and third and fourth times I asked it to do this. Then I went searching for a solution. What I figured out is that in Community Server Dashboard, under Global Settings -> Default Post Settings, there is an option called "Auto Name Posts". If set to "Yes", CS tries to name each post based on the title, rather than use the de facto standard blog post naming convention, which is based on the date the post was created. After I disabled this option, WLW successfully located its post and downloaded my template to my local machine. Success!

My second problem that I've been struggling with is coming up with a quick and easy way to insert formatted, colored code into posts. Before I had WLW, I would cut and paste code and then go and manually add <div> and <span> tags with with embedded CSS styles using the horrible Community Server interface with the TinyMCE editor. Death seemed preferable. If I specified colors like this:

<span style="color: blue;">

things worked just fine. I would get blue text. But if instead I specified:

<span style="color: #0000FF;">

then the TinyMCE editor would modify the tag to look like this:

<span style="color: rgb(0,0,255);">, which would then turn into:

<span>

by the time it was published. I would not get blue text. Argh!

By using WLW with the Code Snippet plugin, as long as I never open up the post in the TinyMCE editor before it's published, everything is fine and my span tags with their color styles are left untouched. Stay away from TinyMCE!

I haven't worked out all of my blog problems just yet. For instance, Code Snippet supports C# code, but not Matlab. But I got two big problems solved, so I will declare success for the moment.

Previously I wrote about choosing ZedGraph as my graphing tool in a C#.Net project for a client. Today I'd like to give an update on how it's going. I have been very pleased with this graphing library, and it has been relatively straightforward to figure out how to make it do what I want it to do.

I'm plotting X-Y data that represents the shape of an airfoil on a 2D plot with grid lines, axes, and axes labels. One of my main concerns was having the ability to force the scales of the X and Y axes to be the same value. If they are locked, my plot will show the airfoil with the right proportions. If not, the airfoil appears to be a lot fatter than it really is.

Figure 1. Airfoil with axis scales locked

Figure 2. Airfoil with axis scales determined independently of one another

The ZedGraph library comes with a ZedGraph control, which contains of a collection of one or more GraphPane objects. The GraphPane holds all of the plotting information for a single plot, including settings and options in addition to the actual plotting data. For my particular problem I only need one GraphPane, which is the default. My first step is to initialize my GraphPane object.

public MainForm()
{
    InitializeComponent();

     // Initialize ZedGraph
    GraphPane graphPane = zedGraphControl1.GraphPane;
    graphPane.XAxis.Title.Text = "Chord";
    graphPane.YAxis.Title.Text = "Thickness";
    graphPane.XAxis.MajorGrid.IsVisible = true;
    graphPane.YAxis.MajorGrid.IsVisible = true;
    graphPane.Legend.IsVisible = false;
    graphPane.Title.IsVisible = false;
    //graphPane.Title.Text = "ZedGraph";
    zedGraphControl1.Resize += new System.EventHandler(this.zedGraphControl1_Resize);    
    graphPane.AxisChangeEvent += new GraphPane.AxisChangeEventHandler(graphPane_AxisChangeEvent);

}

The key here is the last line of initialization, where I add an event handler for the AxisChangeEvent. This event is invoked by calling the AxisChange() method in the ZedGraphControl object. One generally calls this method whenever new plotting data is added to a GraphPane; it takes care of adjusting the minimum and maximum axes values and the grid spacing. AxisChangeEvents are also generated whenever the user zooms or pans a GraphPane. It is not invoked, however, when the user resizes the ZedGraphControl, so I added an event handler for the Resize event that calls the AxisChange() method. The Resize event handler looks like this:

private void zedGraphControl1_Resize(object sender, EventArgs e)
{
    zedGraphControl1.AxisChange();
}

The AxisChange event handler is called after the the ZedGraphControl has recomputed the minimum and maximum axis values and the grid spacing. It adjusts the minimum and maximum values so that the scale in the X and Y directions are equal, and then it recomputes the grid spacing to ensure that the grid does not have too many or too few lines.

private void graphPane_AxisChangeEvent(GraphPane target)
{
    GraphPane graphPane = zedGraphControl1.GraphPane;

    // Correct the scale so that the two axes are 1:1 aspect ratio
    double scalex2 = (graphPane.XAxis.Scale.Max - graphPane.XAxis.Scale.Min) / graphPane.Rect.Width;
    double scaley2 = (graphPane.YAxis.Scale.Max - graphPane.YAxis.Scale.Min) / graphPane.Rect.Height;
    if (scalex2 > scaley2)
    {
        double diff = graphPane.YAxis.Scale.Max - graphPane.YAxis.Scale.Min;
        double new_diff = graphPane.Rect.Height * scalex2;
        graphPane.YAxis.Scale.Min -= (new_diff - diff) / 2.0;
        graphPane.YAxis.Scale.Max += (new_diff - diff) / 2.0;
    }
    else if (scaley2 > scalex2)
    {
        double diff = graphPane.XAxis.Scale.Max - graphPane.XAxis.Scale.Min;
        double new_diff = graphPane.Rect.Width * scaley2;
        graphPane.XAxis.Scale.Min -= (new_diff - diff) / 2.0;
        graphPane.XAxis.Scale.Max += (new_diff - diff) / 2.0;
    }

    // Recompute the grid lines
    float scaleFactor = graphPane.CalcScaleFactor();
    Graphics g = zedGraphControl1.CreateGraphics();
    graphPane.XAxis.Scale.PickScale(graphPane, g, scaleFactor);
    graphPane.YAxis.Scale.PickScale(graphPane, g, scaleFactor);
}

Having the ZedGraph source control was extremely useful in figuring all of this out, particularly for resizing the grid after I have adjust the minimum and maximum axes values. In the end, I did not need to modify any of the ZedGraph source code, but because it is fairly clean and readable, I was able to use it in place of documentation when I needed more explanation than what was provided in the basic tutorial and class descriptions.

So I still feel good about ZedGraph, and my customer is happy, too!

I'm on vacation this week in northern Michigan, enjoying the good weather and time with my family, and also reading quite a bit, including a book called Beautiful Code published by O'Reilly (more about that in another post). One chapter in particular got me fired up to revisit Euler Problem 14, which I started to work on back in May, but got stuck and have been thinking about it ever since.

The problem goes like this:

I'm solving the Euler problems in Matlab, and I'm working on a 867 MHz Mac PowerBook G4. Sometimes I also limit myself to the time it takes to drain the battery, just so that I don't end up spending days and days solving a problem (I have to work once in awhile, and my kids get upset if I ignore them too much :). So my goal is generally to come up with a fast solution, which in most cases rules out brute force.

With these goals in mind, I thought and thought and thought and thought about this problem, trying to come up with a trick that would allow me to rule out most numbers between 1 and 1,000,000. The most obvious thing to try is only the odd numbers. But I wanted to winnow down my search set even further. Not coming up with anything, I gave up and abandoned Euler problems altogether. I had plenty of other things to work on. But this problem kept gnawing at the edges of my brain.

With so many problems in life, sometimes it's best to stop thinking and start doing. So finally this week I threw in the towel and wrote some code to do a brute force search on all odd numbers. It definitely occured to me that once a sequence has been generated for a number, when that number shows up in a subsequent sequence I should reuse what has already been computed. But I got hung up in figuring out how big I should preallocate my array to store the sequence count for each integer, so I even left off that. I just wrote a loop that computes the sequence for every number.

Imagine my surprise when it took only about 15 or 20 minutes to run!

Sometimes brute force isn't so bad. Sometimes the only way to know is to try it.

Here's my code:

function euler14

max_count = 0;
max_nn = 1;
for nn = 3:2:1000000

    new_nn = nn;
    count = 0;

    while new_nn ~= 1
        if mod(new_nn,2) == 0
            new_nn = new_nn/2;
        else
            new_nn = new_nn*3+1;
        end
        count = count + 1;
    end

    if count > max_count
        max_count = count;
        max_nn = nn;
    end

    disp([num2str(nn), ': ', num2str(count)]);
end

disp(['Max: ', num2str(max_count), ' for nn = ', num2str(max_nn)]);

In the past year I've worked on two projects that involved piecing through someone else's Matlab code, making sense of it, making it work, and figuring out how to link to it from a Windows desktop application. In both cases, it was clear that the Matlab had been writen by someone who does not normally spend a lot of time writing code. All of the tenets that we developers hold sacred had been broken: global variables abounded; code was copied over and over; large sections of code were commented out for no apparent reason; and in one of the projects, functions were dispensed with all together in favor of scripts, which often relied on variables that had been declared and set elsewhere in another script. In both cases the original author of the code was not available to explain how things were supposed to work, so fixing the code was a lengthy, tedious process.

It got me wondering...why do engineers write bad Matlab code? Most engineers have had at least one programming course in college, and presumably in that course they are taught how to organize code, declare functions, and pass arguements. Now granted, some people have a harder time than others in figuring out a good architecture for a project. But my observation is that engineers writing Matlab seem to have an aversion to declaring functions. They would rather copy and paste code 20 times than make that code into a function.

And then, as I was working on my most recent Matlab project, it occurred to me what the problem is. When I write C# or C++ code with Microsoft Visual Studio 2005 or 2008, I have IntelliSense, which is a godsend because I can never remember the order in which arguements are supposed to be passed to a function. Matlab does not have anything like IntelliSense, and it has lots and lots of functions, most of them overloaded or with a variable number of arguments. I can't be the only scatter-brained engineer out there who has trouble keeping track of all of my options.

So this blog post is really a plea to The MathWorks folks - give us IntelliSense!

At a minimum, I'd be happy if there were a window that I could have open that would automatically show me the help file for a function as I type its name into an M-file. This would be similar to the context-sensitive help that's available in Labview, another engineering tool that has lots and lots of functions that are hard (for me) to keep track of.

Earlier this month Bill posted his thoughts on software development tools licensing. After several busy weeks of programming for three different clients, I'm finally catching up on reading, but I was too late to post a comment over on Bill's blog, so I'm sharing my thoughts over here.

I'm surprised that the more expensive, less frequently used software development tools don't usually come with floating licenses. The big mechanical engineering packages that I'm familiar with (Nastran, ANSYS, Unigraphics, Catia, Ideas, etc.) can be installed on all users' machines, and then you set up a license server and "check out" a license when you use the tool. If all the licenses are checked out, you bug the guy at the next machine to wrap up his work and quit so you can do your thing. Or you work in the middle of the night, as we all did in graduate school at the end of the semester when everyone was wrapping up projects. At a previous job, when things got particularly busy, we purchased a three month floating license for our finite element software to get us by without us killing one another. In addition, we had the option of checking out the license to a laptop so that we could work when off the company intranet.

Just as Mike Woelmer has experienced with Compuware (see the comments on Bill's post), a single floating license for Matlab or ANSYS or other packages is almost twice as much as getting a single seat, so you have to anticipate several occasional users for it to be worth it. Why haven't the software development tool makers adopted this model, like the mechanical engineering tool makers have?

What is the first triangle number to have more than 500 divisors?

Ok, this one kicked me in the pants, and it kicked my little Mac in the pants as well. I had to relearn how to get divisors quickly and efficiently, which meant seeking the wisdom of others on the Internet. But eventually I got a working program that didn't take all night to run.

For those of you out there doing e-commerce, check out 7 Billion People. They specialize in a more personalized e-commerce experience, and they have a hard-working, experienced team behind them. Take a look at their marketing videos:

A customer asked me to write a very small application that performs some operations on 2D parametric curves. While their code base for their primary software is C++ Win32, I thought this would be a great opportunity to show off C#.Net to them, but that meant coming up with a basic graphics library that would allow me to display (x,y) data as a spline curve or polyline, toggle control points display, and perhaps do some zooming and panning. While the customer is willing to pay for a 3rd party library if the price and feature set are right, I thought I'd look into open source libraries to see what's out there. That's how I discovered ZedGraph.

ZedGraph is billed as a charting library (it can do X-Y, bar and pie charts), but it's X-Y graphing functionality is solid enough to make it appropriate for a basic 2D graphics app. One of my big concerns was the ability to exert control over the scaling along the X and Y axes. I want to be able to force a 1:1 aspect ratio between the X and Y scaling so that the curves that I'm plotting aren't distorted. ZedGraph allows me to set the minimum and maximum values along an axis, but then it takes care of the rest. If I'm showing the axes and a grid behind my curves, ZedGraph automatically figures out where to put tick marks and major and/or minor grid lines.

Naturally I want to be able to control the color and style of the plotted curves, the type of symbol used to represent curve control points, the axes titles, the chart title,  the gridline color, etc. This is very intuitive in ZedGraph, because the architecture of the library and the names of the properties and methods make sense.

Zooming and panning functionality is available by default in the .Net UserControl that comes along with the library. Best of all, the functionality mimics that of many well-known graphics software apps (like ANSYS and Unigraphics).

Finally, the documentation and examples available at the ZedGraph website are decent enough to get the novice beyond setting up a basic graph. I was able to follow the examples to get the look, feel, and functionality my customer required in place within a few hours.

I haven't played around with the library long enough to discover its limitations - I'm sure they're there! But for a very basic 2D graphics app, this library is an excellent choice.

In this problem we are searching for the largest product of four consecutive numbers in any row, column, or diagonal of a given 20 by 20 matrix. This is a nice sort of problem to do in Matlab because it is so well-suited for working with matrices. I wrote a separate function for scooting along a vector looking for the greatest product. Matlab has a built-in diag function that gives you a matrix diagonal as a vector. The 0th diagonal is the central one, the negative indices give you the lower triangle, and positive indices give you the upper triangle. Without further ado...

Sum the primes below 2,000,000.

This is the sort of problem that Matlab handles very well. The code looks pretty, and I got an answer right quick. Now if I type fast, I'll get my code posted before my battery drains and Baseball Tonight ends.