CoCalc provides the best real-time collaborative environment for Jupyter Notebooks, LaTeX documents, and SageMath, scalable from individual use to large groups and classes! Also, H100 GPUs starting at \$2/hour.

Views: 4
Visibility: Unlisted (only visible to those who know the link)
Image: ubuntu2204
Kernel: Python 3 (system-wide)

Python Lanchester's laws simulation

The Lanchester equations are a set of mathematical formulas used to model and analyze the outcome of military engagements between two opposing forces. Developed by Frederick W. Lanchester, a British engineer, during World War I, these equations provide a way to predict the attrition and effectiveness of two battling forces based on their size, firepower, and attrition rates.

The Lanchester equations come in two primary forms: linear and square. The linear form is also known as Lanchester's Law of Primitive Warfare, while the square form is known as Lanchester's Law of Modern Warfare. Here we will only consider the Linear Law (Primitive Warfare). In this model, the effectiveness of each force is assumed to be proportional to its size. It is mostly applicable to ancient or primitive warfare, where engagements were hand-to-hand combats or involved simple projectile weapons like arrows or spears. The linear equations are as follows:

$\frac{dA} {dt} = - \beta B$

$\frac{dB} {dt} = - \alpha A$

The given equations are first-order linear differential equations and can be solved using the method of separation of variables. The equations are of the form:

dx/dt = -ax

Rearrange to isolate terms involving x and t on separate sides of the equation:

dx/x = -a dt

Then integrate both sides with respect to their individual variables:

âˆ« dx/x = âˆ« -a dt

which gives:

ln|x| = -at + C

Solving for x by exponentiating both sides yields:

|x| = e^(-at + C)

Which simplifies to:

x(t) = C e^(-at)

where C is an arbitrary constant that would be determined by an initial condition. The absolute value bars are dropped because C can take any real value, including negative values.

(c) 2022 Marcin "szczyglis" SzczygliÅ„ski

Email: [email protected]

Version: 1.0.0

# R = red force numbers R[t] is the Red force at time = t

# B = blue force numbers B[t] os the Blue force at time = t

# b_s = offensive power of the Blue force

# r_s = offensive power of the Red force

# t = elapsed time

#Here we import Python libraries
import numpy as np                   #import math library
import matplotlib.pyplot as plt      #import plot library
__version__="1.0.0"


The differential equation for the linear case $\frac{dR}{dt}= - \beta \times B$ becomes a difference equation; R(t)-R(t-1) = - $\beta \times B(t-1) \delta t$ or R(t) = R(t-1) - $\beta \times B(t-1) \delta t$

This means the Red force strength, R, at time t-1 is reduced by the 'efficiency factor' $\beta \times$ the Blue force strength B at time t-1. This same equation is used for the reduction in the Blue force.

R.append(Rnew) is a built-in python function that places the value Rnew at the current end of the list R. for t in np.arange(0, T, dt): will loop through the indented statements using the value t that changes during each loop the changes in t are described by the numPy function np.arange(0,T,dt). t will vary from 0 to T in steps of dt. In this case T may be defined as 100 and dt as .t and t would then be, 0, .1 .2 and so on until it reached 100. the loop would then exit to the next statement.

In the next cell we will implement the linear law in python

def calc_linear(R, B, r_l, b_l, t):        #called from linear law
Rnew = R[t] - b_l * B[t]       # compute the next value of the Red force strength
#print("Rnew = ", Rnew, "R[t] = ", R[t])
Bnew = B[t] - r_l * R[t]       # compute the next value of the Blue force strength
# print("Bnew = ", Bnew, "B[t] = ", B[t])
R.append(Rnew)             #put this new value of R into the list R
B.append(Bnew)             #put this new value of B into the list B
#print("calc t = ", t)
return R, B

# linear law
def linear(R0, B0, r_l, b_l, T, dt):

R = [R0,]      #initialize the list R
B = [B0,]      #initialize the list B
for t in np.arange(0, T, dt):
R, B = calc_linear(R, B, r_l, b_l, t)  #call the function calc_linear to update the arrays
if R[t] < 1 or B[t] < 1: break         # continue until one of the forces is eliminated
#    print("t = ", t, "R = ", R[-1], "B =", B[-1])  #print the final values
return R, B

# base parameters:
R0 = 10000  # number of RED units
B0 = 10000  # number of BLUE units

T = 3000  # total number of steps in the simulation
dt = 1     # time interval

r_l = 0.001 # combat efficiency of RED units, each time step each Red kills .1% of Blue
b_l = 0.0012 # combat efficiency of BLUE units, each time step each Blue kills .12% of Red

#Run the linear model.
b_l = 0.001
R, B = linear(R0, B0, r_l, b_l, T, dt) # this line calls the previously defined Linear Law model Python code.

# display result
print("Predicted result of the battleï¼š\n")
if R[-1] > B[-1]:
print("Winner: RED")
else:
print("Winner: BLUE")

# display remaining units info
print("Remaining RED units [", R[-1], "]")
print("Remaining BLUE units [", B[-1], "]")

# display result on plot
t = np.arange(0, len(R) * dt, dt)
#plt.figure(1)
plt.plot(t, R, '--r', label='RED units')
plt.plot(t, B, ':b', label='BLUE units')
plt.xlabel("Time (round)")
plt.ylabel("Number of units")
plt.title("Lanchester's model simulation")
plt.legend()
plt.show()

Predicted result of the battleï¼š Winner: BLUE Remaining RED units [ 497.12393998036214 ] Remaining BLUE units [ 497.12393998036214 ]
print(R[0], R[1], B[0], B[1])

10000 9990.0 10000 9990.0

A number of variations of these equations have been produced and make for interesting reading and exercises. The purpose of this example is to introduce a simple model and subsequent simulation.

What happens if we run this a number of time with a random value of b_1

import matplotlib.pyplot as plt   # import the plotting routines <https://matplotlib.org/>
import random                     # import the routine for generating random numbers <https://docs.python.org/3/library/random.html>
import numpy as np                # numpy is a Python library for scientific computing see <https://numpy.org>
from scipy.stats import gaussian_kde
import seaborn as sns             # another plot library

r_l = 0.001 # combat efficiency of RED units, each time step each Red kills .1% of Blue
b_l = 0.0012 # combat efficiency of BLUE units, each time step each Blue kills .12% of Red


X1 = []   # create a list for b_l
Y1 = []   # create a list for the outcome for B

# run the simulation a number of times
for i in range (0, 10):
b_l = random.random()/1000 + .0005  #create a random number to change b_l randomly
R, B = linear(R0, B0, r_l, b_l, T, dt) # this line calls the previously defined Linear Law model Python code.
X1.append(b_l)    # add the b_l number to its list
Y1.append(B[-1])  # add the final number of B to its list
#   print(b_l)
plt.xlabel("b_l")
plt.ylabel("Number of B units")
plt.title("Lanchester's model simulation results function of b_l")
plt.scatter(X1,Y1)   #plot a scatter diagram of b_l and the resulting values of B
plt.show()
plt.xlabel("b_l")
plt.ylabel("number")
plt.title("Histogram of random number of b_l")
plt.hist(X1)          #plot a histogram of the values of b_l
plt.show()