Battle Simulation With Lanchester Model in C#


This post will show a simple implementation of the famous Lanchester equations to simulate battles “evolution”.

For those that does not know this model, i really recomend this and this PDFs. (especially the first one). For the rest of the post, iam assuming that the reader “WILL read” these recomended texts =P.

Lanchester modeling is pretty beautiful cause it is based in very simple hypothesis and even so, it stills very usefull. (not so good as some stochastic models but good enough for lots of situations ).

In simple ways, we have a system of Differential equations to solve. There are LOOOTS of numeric methods avaliable, i choose one that for me is very intuitive and easy to implement. (using EigenValues and Eigen Vectors -> it is graphically intuitive =P)

The following code shows the Lanchester Model Implementation:
I used plain C# with .net and the alglib for math.
Actually, the alglib was used just to get the eigenvalues and eigenvectors of the lanchester diferential system. It could be easily done by “hand”. (calculate roots of the characteristic equation to find the eigenvalues and then substitute in the formal eigenvector/eigenvalue definition equation to get the eigenvectors =P).

The Squad Class (represents one Army =P)

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace Settlers.Battle
{
    class Squad
    {
        public float numericalStrenght;
        public float effectiveness;
    }
}

The simulator class:

There are some static methods that can be used to discover the winner before the “Battle begin =P” (without explicitelly using the time variable. See the recomended papers to know the reason why we can do this)

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace Settlers.Battle
{
    class LanchesterSimulation
    {
        public LanchesterSimulation(Squad a1, Squad a2)
        {
            double[,] vects = new double[2,2];
            vects[0, 1] = -a2.effectiveness;
            vects[1, 0] = -a1.effectiveness;
            double[] real = new double[2];
            double[] ima = new double[2];
            double[,] left = new double[2, 2];
            double[,] right= new double[2, 2];

            alglib.rmatrixevd(vects, 2, 2, out real, out ima, out right, out left);

            ev1 = real[0];
            ev2 = real[1];

            auv1x = right[1, 0];
            auv1y = right[0, 0];

            auv2x = right[1, 1];
            auv2y = right[0, 1];

            c2 = (a2.numericalStrenght - (a1.numericalStrenght * auv1y) / auv1x ) * auv1x / (-auv2x * auv1y + auv2y * auv1x);
            c1 = (a1.numericalStrenght - auv2x * c2) / auv1x;

        }

        double ev1, ev2;
        double c2, c1;
        double auv1x,auv1y,auv2x,auv2y ;

        public float GetArmy1Size(float t)
        {
            double resp = c1 * auv1x * Math.Exp(ev1 * t) + c2 * auv2x * Math.Exp(ev2 * t);
            return (float)resp;
        }

        public float GetArmy2Size(float t)
        {
            double resp = c1 * auv1y * Math.Exp(ev1 * t) + c2 * auv2y * Math.Exp(ev2 * t);
            return (float)resp;
        }

        public static Squad WhoWillWin(Squad a1, Squad a2)
        {
            float val = a1.effectiveness * a1.numericalStrenght * a1.numericalStrenght - a2.effectiveness * a2.numericalStrenght * a2.numericalStrenght;
            return val > 0 ? a1 : a2;
        }

        public static Squad[] WhoWillWin(Squad[] a1, Squad[] a2)
        {
            float apower = 0;
            float atotal = 0;
            for (int i = 0; i < a1.Count(); i++)
            {
                apower += a1[i].effectiveness * a1[i].numericalStrenght;
                atotal += a1[i].numericalStrenght;
            }

            float bpower = 0;
            float btotal = 0;
            for (int i = 0; i < a2.Count(); i++)
            {
                bpower += a2[i].effectiveness * a2[i].numericalStrenght;
                btotal += a2[i].numericalStrenght;
            }

            float val = apower * atotal - bpower * btotal;
            return val > 0 ? a1 : a2;
        }
    }

}

That is all for today =P