[PYTHON] Calculate the optimal solution to set a world record for decathlon with scipy.optimize

Introduction

Do you know the competition called "decathlon"? In 2 days, 10 events of running (100m, 400m, 1500m, 110mH), jumping (high jump, pole vault, long jump), and throwing (shot put, discus, javelin throw) were performed, and the records of each event were scored and the total record. It is a competition to compete in. The current (June 2020) world record is Kevin MAYER's 9126 points.

By the way, what should I do to get this score efficiently? There seems to be an idea of "scoring evenly in all events", but this is not realistic. This is because there are some events that make it easy to earn points and some that make it difficult to earn points. For example, the record for getting 900 points in each item is as follows. (For reference, Japanese records for each decathlon event are also included)

In order to get 900 points with shot put or discus, you have to break the Japanese record. This is not realistic.

So how is the easiest way to set a world record? Instead of getting only 800 points in shot put, let's get 1000 points in other events, but I don't know if it's easier to get 1000 points in the 100m or 1500m ...

The answer I gave to this question is ** The minimum total score on the "scoring table ** (commonly known as the Hungarian table)" ** "**.

What is a scoring table?

This is a scoring table created by the International Athletics Federation (IAAF), and the calculation formula is different from that of the decathlon. This is made for comparison by event, so you can use it like "The record of 100m 11.00 (886p) that I made is equivalent to 6m83 (886p) in the long jump".

Calculated

We calculated the "minimum score for the 10th event (constraint: the record for the 10th event exceeds 9126 points)".

For example, let's say you aim for a total of 1000 points in the second type of 100m and shot put. The total score when 500 points are given in each event is 972p, The total score at 100m14.00 and shot put 13.32 is 945p, and even if the score is the same, it seems that you can easily earn points. (Actual problem is different ...)

With this kind of feeling, I am looking for how to set a world record with a small score.

Calculation method

The optimization method used minimize in scipy.optimize.

The range of the record was a score of 1p to a world record. (The lower limit of the track event is the world record, and the upper limit of the field event is the world record)

The formula for calculating the score of the decathlon is described in [Wikipedia](https://ja.wikipedia.org/wiki/decathlon #score). I couldn't find the formula for the scoring table officially, but it seems to be calculated by the formula below.

a × (record + b) ^ 2 + c (Coefficients are listed below)

It's surprisingly simple.

The following are the commands used in the calculation.

python

import numpy as np
import pandas as pd
import copy
import sys
from scipy.optimize import minimize, BFGS, LinearConstraint, Bounds
from math import floor
import matplotlib.pyplot as plt

##Coefficient for scoring table
w_score = np.array([
    [24.64221166,-16.99753156,-0.218662048], #100m
    [1.021013043,-78.99469306,0.0029880052], #400m
    [0.0406599253,-384.9950051,0.001205591], #1500m 
    [7.665206128,-25.79302259,0.0141087786], #110mH
    [32.14570816,11.59368894,-5026.080842],  #HJ
    [3.045719921,39.33586031,-4993.213828],  #PV
    [1.931092873,48.34861905,-4993.807793],  #LJ
    [0.0423461436,684.8281542,-19915.72457], #SP
    [0.0040063129,2232.983411,-20003.52492], #DT
    [0.0024031525,2879.797864,-19950.96836]])#JT

##Decathlon coefficient
w_dec = np.array([
    [25.4347,18,1.81],   #100m
    [1.53775,82,1.81],   #400m
    [0.03768,480,1.85],  #1500m
    [5.74352,28.5,1.92], #110mH
    [0.8465,75,1.42],  #HJ
    [0.2797,100,1.35], #PV
    [0.14354,220,1.4], #LJ
    [51.39,1.5,1.05],  #SP
    [12.91,4,1.1],     #DT
    [10.14,7,1.08]])   #JT

##Objective function (scoring table)
def calc_score(x):
    total_score = 0
    for i in range(10):
        total_score += w_score[i,0] * (x[i] + w_score[i,1])**2 + w_score[i,2]
    return total_score

##Constraint function (decathlon score)
def calc_dec(x):
    n = 0
    total_point = 0
    target_point = 9126

    for i in range(10):
        if i in (0,1,2,3):
            total_point += w_dec[i,0] * (w_dec[i,1] - x[i]) ** w_dec[i,2] #100m, 400m, 1500m, 110mH
        elif i in (4,5,6):
            total_point += w_dec[i,0] * (x[i] *100 - w_dec[i,1]) ** w_dec[i,2] #High jump,Long jump,Pole vault
        else:
            total_point += w_dec[i,0] * (x[i] - w_dec[i,1]) ** w_dec[i,2] #Shot put,Discus,Javelin throw
        return_point = total_point - target_point
    return return_point

bounds = Bounds(world_rec[0:4] + min_rec[4:10] , min_rec[0:4] + world_rec[4:10])

cons = (
    {'type': 'ineq', 'fun': calc_dec} 
)

x0 = np.array([10, 46.19, 247.42, 13.59, 1.8, 3.5, 6.06, 10.79, 31.40, 53.67]) #Initial value is appropriate

Almost nothing is described except for the creation of functions. It's very easy to solve an optimization problem with just this ...

result

1. 1. Minimum score when achieving world record (9126 points)

python

result = minimize(calc_score, x0, constraints=cons, method="SLSQP", bounds=bounds)
print(result)
#     fun: 8707.88035324152
#     jac: array([-141.92468262,  -27.99401855,   -5.38891602,  -70.75549316,
#        902.8885498 ,  277.25720215,  221.29797363,   59.95800781,
#         18.4855957 ,   14.31445312])
# message: 'Optimization terminated successfully.'
#    nfev: 569
#     nit: 38
#    njev: 34
#  status: 0
# success: True
#       x: array([ 14.11782497,  65.28575996, 318.7265988 ,  21.17765192,
#         2.45      ,   6.18      ,   8.95      ,  23.12      ,
#        74.08      ,  98.48      ])

100m:14.11 400m:65.28 1500m:318.72 110mH:21.72 High jump: 2.45 (world record) Pole vault: 6.18 (world record) Long jump: 8.95 (world record) Shot put: 23.12 (world record) Discus: 74.08 (world record) Javelin throw: 98.48 (world record)

Total score: 8707

The result was! ww The score increases exponentially, so if you do your best in a field event with a high growth rate Does it mean that even if the track event is appropriate, it will set a world record ...

2. Minimum score when achieving the Japanese record (8308 points)

It's boring if it's just a world record, so I also tried a Japanese record version. The upper limit is "Japanese record for each decathlon".

python

#Only the result is described
result = minimize(calc_score, x0, constraints=cons, method="SLSQP", bounds=bounds)
print(result)
#     fun: 8397.295007256867
#     jac: array([-262.33532715,  -58.31018066,   -7.70007324, -130.22009277,
#        884.24414062,  271.89660645,  216.27709961,   59.32495117,
#         18.29467773,   14.19604492])
# message: 'Optimization terminated successfully.'
#    nfev: 521
#     nit: 30
#    njev: 30
#  status: 0
# success: True
#       x: array([ 11.67464597,  50.4396491 , 290.30574658,  17.29878908,
#         2.16      ,   5.3       ,   7.65      ,  15.65      ,
#        50.23      ,  73.82      ])

100m:11.67 400m:50.43 1500m:290.30 110mH:17.29 High jump: 2.16 (decathlon Japanese record) Pole vault: 5.30 (decathlon Japanese record) Long jump: 7.65 (decathlon Japanese record) Shot put: 15.65 (Decathlon Japan record) Discus: 50.23 (Decathlon Japan record) Javelin throw: 73.82 (Decathlon Japanese record)

Total score: 8397p

The result was that. The tendency to score points in field competitions remains the same ... By the way, the score when Ushiro set a Japanese record is 8609p, so With the above, it will be possible to set a Japanese record efficiently by about 200p w

bonus

Although it has nothing to do with scipy, I tried to graph the score of the decathlon and the calculation formula of the scoring table.

##Graph output
fig, ax = plt.subplots(5, 2, figsize=(14, 25)) 

for i in range(10):
    pltx = np.arange(min(min_rec[i],world_rec[i]), max(min_rec[i],world_rec[i]), 0.01)
    plty = w_score[i,0] * (pltx + w_score[i,1])**2 + w_score[i,2]
    if i in (0,1,2,3):
        plty_p = w_dec[i,0] * (w_dec[i,1] - pltx) ** w_dec[i,2] #100m, 400m, 1500m, 110mH
    elif i in (4,5,6):
        plty_p = w_dec[i,0] * (pltx *100 - w_dec[i,1]) ** w_dec[i,2] #High jump,Long jump,Pole vault
    else:
        plty_p = w_dec[i,0] * (pltx - w_dec[i,1]) ** w_dec[i,2]
    
    ax[i//2, i%2].set_ylim([0,1400])
    ax[i//2, i%2].plot(pltx, plty, color="blue", label="Scoring")
    ax[i//2, i%2].plot(pltx, plty_p,color="orange", label="Decathlon")
    ax[i//2, i%2].set_title(label[i], size=15)
    ax[i//2, i%2].set_xlabel('Record')
    ax[i//2, i%2].set_ylabel('Score / Point')
    ax[i//2, i%2].axvline(x=world_rec[i], ymin=0, ymax=100, ls="--", color="red", label="World Record")
    ax[i//2, i%2].axvline(x=national_dec_rec[i], ymin=0, ymax=100, ls="--", color="green", label="National Decathlon Record")
    ax[i//2, i%2].grid(which = "major", axis = "both", alpha = 0.8,
        linestyle = "--", linewidth = 0.8)
    ax[i//2, i%2].legend(loc='best')

plt.tight_layout()

In the high-scoring area, I feel that the track event has a lower score for the 10 types of scores than the field event (= even if a high score is given, it is difficult to be reflected in the score). I feel like you can understand why we came to the conclusion that "to get a high score with the minimum score, get a high score in the field". ..

At the end

scipy.optimize could be used by me as a statistician amateur. Very convenient ... We hope that this article will increase the number of people (& people who are interested in decathlon) who can easily tackle optimization problems.