Sunday, January 11, 2009

In the last post, I described a fun problem I wanted to play with, and my initial attempt to solve it. I made several ships, figured out how their thrusters worked, and built a simple system for flying. The keys were W for forwards, S for backwards, Q for slide left, E for slide right, A for rotate left, and D for rotate right. Once I flew them around and saw that it was fun, and ship design might be interesting for the player, I wanted to understand better the characteristics of the ships. So I did something that often helps me: visualize the data.

In this diagram you can see the black dots are the possible outputs for the first test ship, viewed from three different perspectives. You can see that it's not completely random. There's a definite shape there. The variety of shapes is what makes this problem interesting. I plotted these for all of my test ships, and learned a lot. For example, look at the lower left plot (with the AD + WS axes). From it we can learn that the ship flies forwards quickly (W) but backwards slowly (S). It can rotate left (A) and right (D) slowly. Going forwards makes it easier to rotate quickly; going backwards makes it harder.

Diagram of random output configurations

While flying the ship I noticed however that there was a second problem I hadn't appreciated at first. When choosing outputs that are “close” to the desired output, a difference of +3 to +4 is nowhere near as significant as a difference of 0 to +0.1. When the player wants zero movement, it's very important to keep it zero. If one key is pressed down, then there should be one non-zero in the output; if two keys are pressed then there should be two non-zeros. The most common case for the random input/output pairs is to have no zeros, but the most common case for the player is to be pressing just one key, and thus needing two zeros!

I often work best when alternating between theory and practice.

Because my random input/output approach was fun but didn't handle the one and two-keystroke cases well, I also looked a few different analytical approaches to compute the reverse mapping directly. I looked at the pseudo-inverse operation, which is like matrix inverse but works with non-square matrices. However, it didn't look like it would help me. I also looked at it as a linear programming optimization problem. That approach looked promising but the Simplex algorithm was more than I wanted to implement.

None of the mathematical approaches I saw directly matched the problem I was trying to solve. It's easy to add the constraint about zeros to a linear programming problem but not to the matrix pseudo-inverse. It might be made to work with the random input/output pairs, by adding some weights to the output components, but the outputs generated from random inputs almost never contained zeros, so I'm never going to get exactly what I want.

One of my habits that seems to work well for me is to alternate between working things out on paper and writing code. The ship thruster physics, the input/output pairs, and the matrix math I first worked on paper, then implemented it. So I went back to paper and pen for this problem. Can I increase the number of outputs that contain some zeros?

First I tried changing the way I picked random inputs, and favored 0 and 1 instead of evenly uniformly choosing any number between 0 and 1. That didn't help much (but I later had an insight related to this change, so maybe it was useful after all). I decided I needed to attack the problem directly. The idea of interpolating between output points was still in my head, and I used that here. If I picked two points on opposite sides of a zero plane, I could find some interpolation that was on the plane. I took pairs of the random input/output pairs and solved for any zeros on the line between them. This worked well! I now had a new set of points that had one zero in the output. By applying the algorithm again, I could find a set of outputs that had two zeros.

In this diagram you can see the black dots are the random outputs projected down to the three planes, and also the white dots, which are formed from finding the point in between two black dots that intersects the plane. The space covered by white dots is more limited than that for black dots. That's because there are some wild movements that can't be controlled if you try to restrict one of the outputs to zero.

Diagram of restricted output configurations

I had up to 1000 random input/output pairs. Solving the equations for every zero plane and every pair that was on opposite sides of the plane meant I needed to solve a matrix equation around 750,000 times. It took a while but it was reasonable. Applying this again to get the outputs with two zeros would've been too slow. And I wanted this to be fast enough so that you could see the flight characteristics as you were editing the spaceship.

Have you ever noticed that you often find something or think of something only after you stop looking for it? I've noticed that this happens for me when solving problems. I had gotten stuck with the brute force technique and I needed something smarter, but I was getting nowhere. So I stopped working on the problem. A few days later, while playing a game, I had an insight that should have been obvious from the start.

All the random inputs are thruster settings from 0.0 to 1.0. We're taking matrix M and multiplying it by a vector T, which must be in a fairly small space. If there are N dimensions, T is an N-dimensional vector, and all its components are between 0.0 and 1.0. That means … the space of all valid thruster configurations is a unit N-dimensional hypercube.

A wise man once told me that it's sometimes easier to solve the general case than the specific case. I had been trying to solve the specific case, of a single input mapped to a single output. And then I'm using computational power to handle lots of them. The general case is to take the entire hypercube and transform it to the output, and see what comes out.

What comes out is a polyhedron in 3 dimensions.

This insight seemed like it came out of nowhere, while I was playing a game, killing orcs. But when I thought about it more, I think it wasn't out of nowhere. If I were smarter I would've figured it out from the very beginning, but some things I had thought touched on the hypercube approach. One clue was that when I picked random numbers for inputs, I found that it useful to bias towards 0.0 and 1.0. My experiments were guiding me towards choosing the vertices of the hypercube. Another clue was that in linear programming, you can think of the valid space as a convex polyhedron in some higher dimensional space. I had spent some time thinking about what the polyhedron represented in my problem but I wasn't able to make the connection. A third clue was that plotting the black and white dots made it look like there were simple polygons involved.

For N thrusters the hypercube has 2N vertices. In this diagram you can see the hypercube vertices, multiplied by the matrix, do cover the space that the random points covered. As before, the black dots are the random sample. The colored dots are the vertices of the hypercube.

Diagram of hypercube vertices

Instead of picking lots of random points, I should pick these points. And then I can take every pair of these and find the interpolation that has a zero output. And then I can take every pair of those and find the outputs with two zeros.

Sometimes I work things out on paper

I worked this out on paper. Page after page of solving equations for test ship 1, following by plotting things out on graph paper, convinced me that this was quite promising.

I was really happy with myself. I had figured out a much more elegant solution that didn't involve random numbers.

I took the equations I solved on paper and wrote some code to solve them. I tried it out on several different ships. And then I ran into a problem. For one of my ships, with 8 thrusters, the algorithm still ran rather slowly. Why? Well, 28 is 256. Finding the one-zero outputs involved examining 215 pairs in 3 dimensions, or around 100,000 systems of equations to solve. It took a few seconds, but I hadn't gotten to the one-zero case yet, which is the really expensive one. For this ship I was almost as slow as the random point approach. The random points were slow but equally slow for all ships, whereas the hypercube vertices got worse every thruster you added. Ick. While starting a cube lamp on my desk, I realized that all the useful zeros should occur along edges, not as interpolations between two arbitrary points. I could reduce 28 pairs X 27 pairs to just 28 pairs X 8. Much faster! However, that was only the one-zero outputs. I still had to run the interpolations again, to find the two-zero outputs. And even with the optimization that was around 1,000,000 systems of equations. I was sure there was another optimization involving surfaces instead of edges but I just couldn't figure out what it was.

So I took another break. I spent some time plotting what I had, to see if there were insights I could draw from the visualization. I plotted sets of points and found that many of them just aren't relevant because they're in the “interior” of the areas. All I care about is the perimeter. To compute the perimeter in two dimensions is a well-known problem: convex hull.

I thought about how to compute a convex hull, and it seemed rather simple. I wrote down my notes, then proceeded to see what the standard algorithms were. None of them seemed as simple as my algorithm, which meant my algorithm probably has a bug. But I couldn't see what it was, so I decided to let it wait a bit. A few days later, still unable to think of a flaw, I sketched out the math on paper, then implemented it. Although I had a few bugs in the code, the algorithm itself worked.

Once I had the convex hull algorithm available, my next step was to use it to simplify the one-zero outputs into a convex hull before finding the two-zero points.

In this diagram you can see the points that make up the convex hull. As before, the white dots are the interpolations between random points, to find positions on the zero plane. The colored points with black borders around them are the interpolations between hypercube vertices, fed into the convex hull algorithm. You can see that all the points I computed from the random set lie within the convex hull computed from the hypercube vertices.

Diagram of convex hulls

In fact the two-zero points always occur between adjacent points in the convex hull, so I had gotten yet another optimization out of this change. However, I discovered that my convex hull algorithm does have a flaw. When there are collinear points it doesn't always remove all of them. The flaw remains in my code; it doesn't seem to be a big deal in practice, because my code is now “fast enough”. For each of the five test ships, generating the flight parameters takes less than a second.

Another thing that wasn't clear to me when I started was that when you press both W and Q, you want to go both left and forwards, but there's not a fixed ratio between the two. Sometimes you can get both and sometimes you have to make a tradeoff, and none of the techniques worked out of the box for that. I think this might be something to leave to the ship designer, once the flight characteristics are known. For example, in the above diagrams you could have a way to mark what you want a combination of keys to do.

It's been a fun journey. I've taken off some rust from my skills. And I learned a bit more about linear programming, the simplex algorithm, physics, and convex hulls.

However, I haven't yet written an editor. But I now have the physics and math worked out, so once I have an editor, you'll be able to see the flight characteristics of the ship you're editing. I'll post the demo and source here. It will be nice if there are interesting tradeoffs for the player. For example, if you put jets on the wingtips, you can rotate very quickly, but perhaps it adds a lot of mass to make the wings sturdy enough to take jets, and that extra mass slows you down in terms of forward acceleration. Once I have the editor I'll get a better feel for whether there are interesting tradeoffs to be made, and whether this would be useful to have in a game.

Labels: , ,

9 comments:

Matthew Skala wrote at January 11, 2009 5:25 PM

It seems to me that you're doing an awful lot of computational work to avoid a little bit of mathematics. The basic problem is that you have a matrix M representing the effects of the thrusters, and a vector t representing which thrusters you're going to activate, and you want to find a value of t such that M*t is of a certain form - usually, so that it's axis-aligned (has only one nonzero entry). Is that it?

Let's take this in steps. Suppose you had exactly as many thrusters as entries in your vector, so M is square, and you could activate any thruster in any amount, including negative. Then you just invert the matrix and you're *done*. You get a new matrix that contains a vector to produce each thrust pattern of the form [1 0 0], [0 1 0], ... It is possible for some matrices to not produce any answer, if your thrusters are in bad places; those would happen to be matrices of zero determinant, but it doesn't really matter what that means; you can just detect it and tell the designer "your design sucks!"

Next step: suppose you have more thrusters, for instance four thrusters and only three dimensions of output. That's kind of nice because it reduces the chance that your thrusters will turn out to be in bad places (you can probably use the extra thruster to break the "tie"); but it means that for any given output there may be multiple thruster activation patterns that all produce that output, so you don't get a unique solution. What people normally do in that kind of case is least-squares: they say "Because I have to choose among many different solutions, I will prefer the one that has the least sum of squares of the variables of my vectors." In your case that's much like saying "Among different thruster activation patterns that achieve the result, I will prefer the solution that minimizes the total amount of thruster fuel spent" which seems like a reasonable way of doing it. That's also easy to do mathematically - none of this Monte Carlo baloney. What you need to do is find the Moore-Penrose matrix inverse, which is one of the pseudo-inverses you were looking at. You can find a formula for computing that in Mathworld; note that the superscript H means "conjugate transpose", which is just transpose because you're not using complex numbers. Feed in your thruster matrix M, you get out an inverse matrix P, with the property that if M*t=x (say x=[0 1 0] or whatever), then P*x is the smallest value of t for which that holds, where "smallest" means "the sum of the squares of the elements is minimized".

Presumably the thrusters have a maximum activation level they can handle. The easy thing to do in that case is take one of the previous solutions, and then scale the result so that the most-activated thruster is activated to its maximum. Then you also get a nice result telling you just how much thrust you can get out of the system in each dimension of thrust - and the designer would probably like to know about that.

But what if you can't activate your thrusters to a negative degree, and what if they all have different limits on how much they can be activated and you don't want to fake that by scaling the input? Then I think you need to set up a linear programming problem for each dimension. It'll be along the lines of "Minimize total thruster activation (which could be weighted if some are more expensive than others) subject to M*t=[0 1 0], and no thruster activated less than zero, and no thruster activated more than its maximum." You can solve that with the simplex method, but DON'T implement it yourself. Go download a library or something. If you must implement it yourself, go get a copy of Numerical Recipes and follow the recipe; you shouldn't have to understand it to implement it, and you shouldn't try to reinvent it.

There are even better ways to solve linear programming problems (interior-point methods), but their advantages don't show up for the problem size you're looking at and so they aren't worth it.

I think almost any solution that requires you to test a large number of guesses at random or educated-random, is almost certainly the Wrong Answer. That's heavy computational artillery, suitable for cases where (for instance) there's a complicated nonlinear relationship among the different thrusters so that activating A and B is different from the sum of activating either alone. I don't think that heavy artillery can be justified for a simple linear problem like this one.

Note, also, that your "convex hull" solution is going to *be* the simplex method by the time you're done tinkering with it. How do you think the simplex method works? It pretty much exactly comes down to building a polytope representing the constraints (which you're calling the convex hull of the projected hypercube) and then exploring its vertices and surfaces to find the best solution. You're taking a long and unnecessarily complicated route to get to that.

Amit wrote at January 11, 2009 6:03 PM

Thanks Matthew! I hadn't made the connection to least squares but that makes a lot of sense. My math is quite rusty and least squares hadn't come to mind.

Monte Carlo is cheap (~15 lines of code in this case), and it's quite useful for getting a quick sense of what the data looks like, especially when it may not be linear. I started with it and then looked for other approaches. Linear programming seemed like the best fit among the things I looked at. However, I didn't find any Simplex algorithm implementations in Actionscript, and I didn't want to implement it myself, which is why I went looking for something else. I disagree that convex hull will become simplex. Finding a convex hull in 2 dimensions seem quite a bit simpler than the simplex algorithm.

Whether you go for the computational approach or the analytical approach depends on how you value computer time and your time. I tend to start with approaches that make the computer do lots of work, and I then look for other approaches if it doesn't work well. My goal here isn't to solve the math or find the optimal solutions; it's to write a quick and dirty spaceship editor demo, to see whether it'd be interesting for a game. If Monte Carlo was satisfactory, I would've stuck with it. I had hoped to finish this project in a few days but it took a bit longer than that, because Monte Carlo didn't work out.

Fusun wrote at January 12, 2009 6:04 AM

Hello

You might be interested in http://fusun.ch/cms/index.php?id=78.

greets

Amit wrote at January 12, 2009 10:11 AM

Fusun: looks interesting! Especially the documents page.

Amit wrote at January 12, 2009 12:30 PM

I took another look at pseudo-inverses and I think I can't use them for this problem. For example, in test ship 1, the “best” way to slide left is thruster input [ 0.0625 0.1875 1 0 ]. But the pseudo-inverse tells me that I need [ -0.0625 0.0625 0.5 -0.5 ]. I'm not sure how to impose the non-negative thruster value constraint when using pseudo-inverses. Linear programming seems like a better fit.

Roundedge wrote at January 19, 2009 10:39 PM

I've been considering this exact same concept, however I hadn't even started to chew on the specifics. I've been more interested in what implications it could have for gameplay. One thing that interests me is that the structure of the ship can be encoded as a genotype, and then you can apply genetic algorithms to create ships which meet certain parameters. For instance, you could have a game, maybe similar to geometry wars, which could slowly evolve better ships suited to your playing style. Or you could have a space RPG which spans many generations, and the ships in the galaxy evolve as the technology gets better.

Darren R wrote at January 20, 2009 9:10 PM

Hi Amit,
I have spent some time considering this problem since I read your blog yesterday. The idea seems to have great gaming potential. Considering Matthew Skala's comments, I disagree (as you have shown by example) that the pseudo-inverse is appropriate. There is just no way to add constraints such as 'no negative thrust' into linear algebra routines. I think that the simplex method is the way to go to provide quick mappings for your controller.
The key is to shape the constraints correctly. If we follow Wikipedia then this problem should probably take the form of a 'dual' problem, that is:

minimise THRUST
subject to M*T = F
and T >= 0

The key here is to recognise that the strict equality condition can be broken up into two parts:
M*T >= F
-(M*T <= F) i.e. -M*T >= -F

By including each of these constraints in the simplex tableau I think you will get the answer.

Of course this will take some time to implement unless you can find a simplex code on the web

Amit wrote at January 23, 2009 8:56 AM

Hi Darren,

Yes, linear programming is the best match for this problem; I just didn't want to implement Simplex myself. The common cases are 2d planes and so I used convex hulls instead of Simplex for those, and then I use an approximation for the uncommon case (all 3 dimensions).

I think if you're using this for a real game, you'd probably want to use Simplex. For the quick prototype I am working on, the goal is to figure out whether this is a fun/interesting idea, and the approximation doesn't change that. I'm not even convinced it was worth spending the effort to find the intersections with zero planes, and convex hulls. I should do that after I figure out whether it's fun, and only if it turns to be fun. :)

Matthew Skala wrote at January 28, 2009 6:50 PM

Do note that I said matrix pseudo-inverse *if* you were okay with activating a thruster to a negative degree, and that you probably wouldn't be. I don't want to be blamed for a flat "Use matrix pseudo-inverse, it's what you want" statement.

One reason I like using a ready-made linear programming solver is that then you can add extra constraints to model other things as your problem becomes more complicated. That's quite hard to do if you're building a custom "like linear programming, but not really admitting that that's what it is" solver.