Exemple d'application de Python à l'ingénierie: recherche de paramètres de formule de calcul approximative à partir de données réelles

Lorsque vous souhaitez mesurer la température, etc. avec un appareil intégré, il existe des cas où vous ne pouvez pas réellement mesurer à la position souhaitée et vous devez la calculer et l'approcher d'une manière ou d'une autre. puis,

――Pouvez-vous faire une bonne formule d'approximation? --Comment déterminer les paramètres de la formule

C'est un problème.

Dans cet article, je vais vous montrer un exemple de la façon de déterminer facilement ce dernier paramètre en utilisant scipy en Python.

――Enfin, il est censé être appliqué aux périphériques embarqués qui ne peuvent effectuer qu'une mauvaise arithmétique entière. --Traitez le système que vous souhaitez approximer comme une boîte noire dans une certaine mesure.

Toutes les opérations suivantes sont effectuées sur le bloc-notes Jupyter. Il est publié sur github.

• daisukelab/general-learning/Python example - simple parameter optimization.ipynb

Lire les données

En supposant que les données ci-dessous sont prêtes sous forme de fichier CSV, commencez par lire ces données.

--t fois --x Valeurs mesurées telles que les capteurs réellement obtenus --y La valeur que vous voulez réellement obtenir, vous voulez approximer cette valeur à partir de x

Le paramètre de problème est "y ne peut pas être mesuré en fonctionnement réel, il est donc approximé à partir de la valeur de x". (Bien que y puisse être mesuré par une méthode spéciale sur le site de développement, il arrive souvent que seul x puisse être obtenu au stade de l'exploitation.)

Correspond aux données de l'enseignant d'apprentissage automatique y et aux données d'entrée x.

%matplotlib inline
import matplotlib.pyplot as plt
import pandas
from scipy import optimize
import numpy as np

#lire csv
datafile = 'xy_data.csv' #Nom du fichier de données
df = pandas.read_csv(datafile)
#Visualisons
df.plot()
plt.show()
#J'afficherai également la valeur numérique des données(100 à 5)
df[100:105]

x(オレンジ)を使ってy(緑)を近似したい。

Utilisons maintenant x pour approcher y.

Modèle 1-Approximation avec équation linéaire

y peut être vu comme 1.? Temps de x. Par conséquent, nous avons décidé d'approximer avec f (x) = a * x + b et de trouver les paramètres a et b avec scipy.

#Définir une approximation
def f1(x, a, b): #Lister les paramètres requis après une entrée
    return a*x + b

#Résumer les calculs approximatifs en fonctions
def param_solver(f, x, y, debug_print=True):
    '''Fonction f(x)Calculez les paramètres ajustés avec la valeur initiale supposée de sorte que l'entrée x corresponde à y.'''
    params, params_covariance = optimize.curve_fit(f, x, y) # <----------Calcul d'optimisation ici
    if debug_print:
        print('params =', params) #Paramètres obtenus
        print('params_covariance =', params_covariance) #Covariance avec ce paramètre
        print('standad deviation errors = ', np.sqrt(np.diag(params_covariance))) #écart-type
    return params

#Convertir les données en données numpy une fois
x = np.array(df['x'].tolist())
y = np.array(df['y'].tolist())
t = np.array(df['t'].tolist())

#Calculer les paramètres
params = param_solver(f=f1, x=x, y=y) 

params = [ 1.24075239e+00 2.31148413e+03] params_covariance = [[ 3.47128802e-04 -4.06799576e+00] [ -4.06799576e+00 5.46259540e+04]] standad deviation errors = [ 1.86313929e-02 2.33721959e+02]

Vous pouvez facilement trouver les paramètres approximatifs dans *** ʻoptimize.curve_fit` *** dans scipy.

Visualisez les approximations avec les paramètres trouvés

À quel point cela pourrait-il être?

#Combinez les visualisations en fonctions
def plot_all(f, x, y, params):
    fig, ax = plt.subplots()
    ax.plot(t, f(x, *params), 'b-', label="f(x)")
    ax.plot(t, x if len(np.array(x).shape) == 1 else x[0], 'r-', label="x")
    ax.plot(t, y, 'g-', label="y")
    ax.legend()
    plt.show()

plot_all(f1, x, y, params)

C'était impossible car je viens de multiplier x (rouge) par une constante par rapport à y (vert) ... Changeons le modèle.

Équation quadratique introduite par le modèle 2

Approximative avec f (x) = a * x ^ 2 + b * x + c.

#Définir une approximation
def f2(x, a, b, c):
    return a*x**2 + b*x + c

#Calculer les paramètres
params2 = param_solver(f=f2, x=x, y=y) #Faites deviner trois éléments

#Visualiser
plot_all(f2, x, y, params2)

params = [ -5.82931117e-05 2.42030720e+00 -2.33839533e+03] params_covariance = [[ 1.78832431e-11 -3.61865487e-07 1.42649657e-03] [ -3.61865487e-07 7.56116875e-03 -3.16641976e+01] [ 1.42649657e-03 -3.16641976e+01 1.51375869e+05]] standad deviation errors = [ 4.22885837e-06 8.69549812e-02 3.89070519e+02]

C'est assez similaire, mais les sommets ne sont pas très proches. Si vous regardez de plus près, il semble que les points de la courbe sont différents dans le temps entre x et y. Considérez également le facteur temps.

Modèle 3-Introduire la différenciation

Approximative en utilisant le terme différentiel dx / dt avec f (x) = a * x ^ 2 + b * x + c * dx / dt + d

def f3(xs, a, b, c, d):
    x, xdelayed = xs
    dx = x - xdelayed #La différenciation est exprimée par la différence avec la composante de retard
    return a*x**2 + b*x + c*dx + d;

def make_delayed(d, delay_step):
    '''Premier retard_étapes d[0]Et le retard arrière_Créer un d avec un pas tronqué'''
    d = list(d)
    n = len(d)
    result = np.array([d[0] for i in range(delay_step)] + list(d[:-delay_step]))
    return result

#Créer des données avec 10 heures de retard
x10 = make_delayed(x, delay_step=10)

#Calculer et visualiser les paramètres
params3 = param_solver(f=f3, x=(x, x10), y=y) #Donner x deux entrées
plot_all(f3, (x, x10), y, params3)

params = [ -3.54499345e-05 2.03961022e+00 3.95514919e+00 -2.84616185e+03] params_covariance = [[ 2.27985415e-12 -4.55332708e-08 2.90736635e-08 1.64730881e-04] [ -4.55332708e-08 9.39580914e-04 -4.84532248e-04 -3.67720697e+00] [ 2.90736635e-08 -4.84532248e-04 5.03391759e-03 -6.46259873e-01] [ 1.64730881e-04 -3.67720697e+00 -6.46259873e-01 1.79598365e+04]] standad deviation errors = [ 1.50991859e-06 3.06525841e-02 7.09501063e-02 1.34014314e+02]

Cependant, la f (x) approximée est légèrement différente de y (vert) dans la seconde moitié. Ici, j'ai utilisé un simple x (t) --x (t-10) comme méthode de différenciation, mais ce délai de 10 temps était-il bon?

Je rechercherai également cette composante de retard (ci-après td) en tant que paramètre.

Trouvez la valeur optimale pour l'approximation différentielle

#Une fonction d'évaluation qui quantifie la qualité de la temporisation td avec la cible et f(x)Valeur absolue de la différence de(Distance L1)Défini comme la somme de
def cost(td, *args):
    _f, _x, _y = args
    #Créer des données de retard
    xdelay = make_delayed(d=_x, delay_step=td)
    #Calculer l'optimisation
    _params = param_solver(f=_f, x=(_x, xdelay), y=_y, debug_print=False)
    #Calculer la somme des distances
    return np.sum(np.abs(_y - _f((_x, xdelay), *_params)))

print('cost(td=5) = ', cost(5, f3, x, y))
print('cost(td=10) = ', cost(10, f3, x, y))

cost(td=5) = 178246.4667 cost(td=10) = 149393.276741

Tout d'abord, j'ai décidé de la fonction d'évaluation facilement, et lorsque j'ai pris la valeur avec td = 5, 10 comme essai, j'ai pu confirmer quantitativement que 10 est meilleur que 5.

Encore une fois, scipy optimise.

def cost_for_brute(td, *args):
    '''comportement de la brute(Enfin td est donné dans la liste)Wrapper selon, toujours donner une valeur entière au coût'''
    try:
        td = int(td[0])
    except:
        pass
    return cost(td, *args)

#chercher
result = optimize.brute(cost_for_brute, [slice(1, 50, 1)], args=(f3, x, y), full_output=True)
final_td, final_cost, tds, costs = result #Td optimal,Coût optimal,Plage de recherche td,Chaque coût
final_td = int(final_td[0]) # [23.]Puisqu'il est retourné sous forme de liste de nombres réels comme, il est converti en entier.
print('final optimum td = ', final_td, 'for cost = ', final_cost)

#Visualisation
fig, ax = plt.subplots()
ax.plot(tds, costs, 'b-', label="cost(td)")
ax.plot(final_td, final_cost, 'go', label="Minima")
ax.legend()
plt.show()

final optimum td = 23 for cost = 123800.706755

Vous pouvez voir qu'il n'y a qu'une seule solution optimale qui est bien convexe et qui a la plus petite valeur de fonction d'évaluation. Maintenant, visualisons l'état d'approximation avec ce paramètre td.

def solve_and_plot_by_td(td, f, x, y, debug_print=True):
    #Créer des données de retard
    xdelay = make_delayed(d=x, delay_step=td)
    #Calculer l'optimisation
    params = param_solver(f=f, x=(x, xdelay), y=y, debug_print=debug_print)
    #Visualisation
    plot_all(f, (x, xdelay), y, params)
    return params

f3_params = solve_and_plot_by_td(final_td, f3, x, y)

params = [ -3.56938479e-05 2.02134829e+00 1.78433928e+00 -2.62914982e+03] params_covariance = [[ 1.41335136e-12 -2.83473596e-08 7.69752574e-09 1.03708164e-04] [ -2.83473596e-08 5.86738635e-04 -1.35889230e-04 -2.30772756e+00] [ 7.69752574e-09 -1.35889230e-04 6.07763010e-04 -9.90337036e-02] [ 1.03708164e-04 -2.30772756e+00 -9.90337036e-02 1.11544635e+04]] standad deviation errors = [ 1.18884455e-06 2.42226884e-02 2.46528499e-02 1.05614693e+02]

Il s'est amélioré. La valeur de la fonction d'évaluation est spécifiquement optimisée comme suit, et f (x) chevauche presque la cible y sur le graphique.

cost(td=10) =  149393.276741
cost(td=23) =  123800.706755

Si l'objectif était atteint avec cela, l'expression approximative f3 () était-elle vraiment optimale?

Par exemple, le terme quadratique «a * x ** 2» était-il nécessaire? Il y a aussi une différenciation, mais un terme d'intégration n'est-il pas nécessaire?

Modèle 4-Omettre le terme quadratique

Approximative avec f (x) = b * x + c * dx / dt + d, en omettant x ^ 2.

def f4(xs, b, c, d):
    x, xdelayed = xs
    dx = x - xdelayed #La différenciation est exprimée par la différence avec la composante de retard
    return b*x + c*dx + d;

#Implémentation de la recherche
result = optimize.brute(cost_for_brute, [slice(1, 50, 1)], args=(f4, x, y), full_output=True)
final_td, final_cost, tds, costs = result #Td optimal,Coût optimal,Plage de recherche td,Chaque coût
final_td = int(final_td[0]) # [23.]Puisqu'il est retourné sous forme de liste de nombres réels comme, il est converti en entier.
print('final optimum td = ', final_td, 'for cost = ', final_cost)

#Visualisation
fig, ax = plt.subplots()
ax.plot(tds, costs, 'b-', label="cost(td)")
ax.plot(final_td, final_cost, 'go', label="Minima")
ax.legend()
plt.show()

final optimum td = 32 for cost = 251326.900162

J'ai trouvé que le coût était mauvais. Au cas où, visualisons l'état d'approximation.

solve_and_plot_by_td(final_td, f4, x, y)

params = [ 1.28868602 1.44297502 176.91492786] params_covariance = [[ 6.30546320e-05 3.59497597e-05 -7.78121379e-01] [ 3.59497597e-05 1.08222102e-03 -1.60091104e+00] [ -7.78121379e-01 -1.60091104e+00 1.21028825e+04]] standad deviation errors = [ 7.94069468e-03 3.28971279e-02 1.10013101e+02]

Au premier coup d'œil, vous pouvez voir qu'il s'est détérioré dans le graphique. Comme mentionné ci-dessus, je comprends que le terme quadratique était nécessaire.

Ensuite, introduisons le terme d'intégration.

Modèle 5-Introduire le terme d'intégration

Approximative en utilisant le terme d'intégration «somme (x)» avec «f (x) = a * x ^ 2 + b * x + c * dx / dt + d * somme (x) + e».

#Obtenez la valeur intégrée dans la section mobile. Divisez par le nombre de données à ajouter pour éviter le dépassement de chiffres=>∴ Moyenne mobile ...
def integral(d, span):
    d = list(d)
    n = len(d)
    dspan = [d[0] for i in range(span - 1)] + d #Span au début-Préparez un tableau avec un ajouté
    return np.array([sum(dspan[i:i+span])/span for i in range(n)])#Créez des données en ajoutant des intervalles lors du déplacement

#Définir une approximation
def f5(xs, a, b, c, d, e):
    x, xdelay, xsum = xs
    dx = x - xdelay
    return a*x**2 + b*x + c*dx + d*xsum + e

#Une fonction d'évaluation qui quantifie la qualité du temps d'intégration ti, la cible et f(x)Valeur absolue de la différence de(Distance L1)Défini comme la somme de
def cost5(ti, *args):
    f, x, xdelay, y = args
    #Créer des données intégrées
    xsum = integral(x, ti)
    #Calculer l'optimisation
    params = param_solver(f=f, x=(x, xdelay, xsum), y=y, debug_print=False)
    #Calculer la somme des distances
    return np.sum(np.abs(y - f((x, xdelay, xsum), *params)))

#Créer des données de retard(* Utilisez la valeur optimale obtenue la dernière fois.)
xdelay = make_delayed(d=x, delay_step=final_td)
#Essayez ti=Calculez le coût lorsque 10
print(cost5(5, f5, x, xdelay, y))
print(cost5(10, f5, x, xdelay, y))

105806.884719 105436.131801

def cost5_for_brute(td, *args):
    '''comportement de la brute(Enfin td est donné dans la liste)Wrapper selon, toujours donner une valeur entière au coût'''
    try:
        td = int(td[0])
    except:
        pass
    return cost5(td, *args)

#Implémentation de la recherche
result = optimize.brute(cost5_for_brute, [slice(1, 400, 1)], args=(f5, x, xdelay, y), full_output=True)
final_ti, final_cost, tis, costs = result #Ti optimal,Coût optimal,Plage de recherche ti,Chaque coût
final_ti = int(final_ti[0]) # [23.]Puisqu'il est retourné sous forme de liste de nombres réels comme, il est converti en entier.
print('final optimum ti = ', final_ti, 'for cost = ', final_cost)

#Visualisation
fig, ax = plt.subplots()
ax.plot(tis, costs, 'b-', label="cost5(ti)")
ax.plot(final_ti, final_cost, 'go', label="Minima")
ax.legend()
plt.show()

final optimum ti = 47 for cost = 89564.819536

Vous pouvez voir que la fonction d'évaluation a cette fois trois solutions locales de l'ordre de [1 400). Il s'est avéré que c'était la solution optimale 47.

def solve_and_plot_by_ti(ti, f, x, xdelay, y, debug_print=True):
    #Créer des données intégrées
    xsum = integral(x, ti)
    #Calculer l'optimisation
    params = param_solver(f=f, x=(x, xdelay, xsum), y=y, debug_print=debug_print)
    #Visualisation
    plot_all(f, (x, xdelay, xsum), y, params)
    return params

f5_params = solve_and_plot_by_ti(final_ti, f5, x, xdelay, y)

params = [ -3.03782209e-05 9.91815249e+00 -4.26387796e+00 -7.99927214e+00 -2.47884462e+03] params_covariance = [[ 6.63718559e-13 3.34745759e-08 -2.98614480e-08 -4.66996212e-08 4.64742006e-05] [ 3.34745759e-08 7.18603686e-02 -5.03435446e-02 -7.23475214e-02 -1.40511200e+00] [ -2.98614480e-08 -5.03435446e-02 3.54804571e-02 5.08234965e-02 2.46752985e-01] [ -4.66996212e-08 -7.23475214e-02 5.08234965e-02 7.31066200e-02 3.71334580e-01] [ 4.64742006e-05 -1.40511200e+00 2.46752985e-01 3.71334580e-01 4.98091129e+03]] standad deviation errors = [ 8.14689241e-07 2.68067843e-01 1.88362568e-01 2.70382359e-01 7.05755715e+01]

La somme des distances L1 renvoyées par la fonction d'évaluation semble toujours aussi grande que 89564, mais le graphique a atteint le point où l'expression approximative correspond presque à la cible. En introduisant le terme d'intégration, c'est la récolte que le délai d'approximation est amélioré.

Formule d'approximation optimale jusqu'à présent

def f5(xs, a, b, c, d, e):
    x, xdelay, xsum = xs
    dx = x - xdelay
    return a*x**2 + b*x + c*dx + d*xsum + e

D'autre part, il a été constaté que chaque paramètre doit utiliser les valeurs suivantes.

params = [ -3.03782209e-05   9.91815249e+00  -4.26387796e+00  -7.99927214e+00 -2.47884462e+03]

Enfin, considérons deux autres.

• Avez-vous vraiment besoin du terme x ^ 2 ??
• Lorsqu'il est appliqué à des périphériques intégrés, il existe des cas où 50 tampons de moyenne mobile ne peuvent pas être conservés. Que faire dans ce cas?

Modèle 6 - Réaffirmer le besoin de termes secondaires

En regardant la valeur du paramètre a = params [0], c'est "-3.03782209e-05".

En d'autres termes, le terme quadratique *** contribue peu au résultat et n'est pas nécessaire ?? ***. Modifions réellement la formule et vérifions-la.

#Définir une approximation
def f6(xs, b, c, d, e):
    x, xdelay, xsum = xs
    dx = x - xdelay
    return b*x + c*dx + d*xsum + e

#chercher
result = optimize.brute(cost5_for_brute, [slice(1, 400, 1)], args=(f6, x, xdelay, y), full_output=True)
final_ti, final_cost, tis, costs = result #Ti optimal,Coût optimal,Plage de recherche ti,Chaque coût
final_ti = int(final_ti[0]) # [23.]Puisqu'il est retourné sous forme de liste de nombres réels comme, il est converti en entier.
print('final optimum ti = ', final_ti, 'for cost = ', final_cost)

#Visualisation
fig, ax = plt.subplots()
ax.plot(tis, costs, 'b-', label="cost(ti)")
ax.plot(final_ti, final_cost, 'go', label="Minima")
ax.legend()
plt.show()

solve_and_plot_by_ti(final_ti, f6, x, xdelay, y)

final optimum ti = 339 for cost = 152778.51452

params = [ 1.45958277 1.24676966 -0.27112565 246.80198745] params_covariance = [[ 1.74462598e-04 -1.29255162e-04 -2.11185881e-04 -4.55806349e-01] [ -1.29255162e-04 8.85113218e-04 2.42461052e-04 -1.11227411e+00] [ -2.11185881e-04 2.42461052e-04 3.35043878e-04 -8.63631214e-02] [ -4.55806349e-01 -1.11227411e+00 -8.63631214e-02 7.95856480e+03]] standad deviation errors = [ 1.32084291e-02 2.97508524e-02 1.83042038e-02 8.92107886e+01]

*** Après tout, le terme quadratique était nécessaire. *** ***

Si la valeur du paramètre est petite, vous pouvez penser que la contribution est faible, mais il est clair que vous ne pouvez pas porter de jugement tant que vous ne l'avez pas vérifié.

Enfin, pour un petit système embarqué qui ne peut pas prendre une grande valeur pour ti (ne peut pas prendre une grande mémoire tampon d'intégration, ne peut pas gérer de grands entiers, est difficile à calculer avec plusieurs bits), regardez ce qui se passe lorsque ti est petit. regarder.

Lors de l'utilisation du modèle 5 avec un petit ti

Il existe jusqu'à 10 solutions locales pour ti. Et cette affaire?

#Implémentation de la recherche
result = optimize.brute(cost5_for_brute, [slice(1, 30, 1)], args=(f5, x, xdelay, y), full_output=True)
final_ti, final_cost, tis, costs = result #Ti optimal,Coût optimal,Plage de recherche ti,Chaque coût
final_ti = int(final_ti[0]) # [23.]Puisqu'il est retourné sous forme de liste de nombres réels comme, il est converti en entier.
print('final optimum ti = ', final_ti, 'for cost = ', final_cost)

#Visualisation
fig, ax = plt.subplots()
ax.plot(tis, costs, 'b-', label="cost(ti)")
ax.plot(final_ti, final_cost, 'go', label="Minima")
ax.legend()
plt.show()

solve_and_plot_by_ti(final_ti, f5, x, xdelay, y)

final optimum ti = 9 for cost = 105290.231346

params = [ -3.35077051e-05 6.34533112e+00 8.26308638e-01 -4.35919543e+00 -2.65976841e+03] params_covariance = [[ 9.91381369e-13 8.29449613e-10 1.80437951e-09 -2.06410309e-08 7.13151940e-05] [ 8.29449613e-10 4.55158040e-02 -4.96419025e-03 -4.52708548e-02 -3.90835831e+00] [ 1.80437951e-09 -4.96419025e-03 7.59625712e-04 4.90797861e-03 2.31789824e-01] [ -2.06410309e-08 -4.52708548e-02 4.90797861e-03 4.54355845e-02 2.30925244e+00] [ 7.13151940e-05 -3.90835831e+00 2.31789824e-01 2.30925244e+00 7.83076217e+03]] standad deviation errors = [ 9.95681359e-07 2.13344332e-01 2.75613082e-02 2.13156244e-01 8.84915938e+01]

Elle n'est pas si inférieure à la solution optimale lorsque «ti = 47».

Lorsque la mise en œuvre est difficile, opérer avec «ti = 9» est souvent un compromis.

Résumé

Après tout, j'ai pu faire une bonne approximation avec PID (proportionnel / intégral / différentiel) + terme quadratique.

Étant donné que la cible initiale de cette époque était un sujet similaire au système de contrôle classique, je pense que c'est une réponse attendue pour ceux qui vont dans cette direction. Lorsque vous pouvez penser à une telle formule d'approximation, si vous utilisez bien scipy en Python, vous pouvez facilement vérifier et trouver un paramètre approprié.

Aussi, à ce moment-là, je pense qu'il est le plus rapide de laisser des notes techniques avec des données et des formules de calcul à l'aide de Jupyter Notebook.

Supplément et développement futur

Cette fois, les données ont été expliquées dans un échantillon. Cela peut être bien si le système ne fluctue pas autant, mais en réalité, il est plus sûr de procéder statistiquement en prélevant 30 ~ échantillons.

À ce moment-là, il y a encore du travail à faire, par exemple s'il faut faire la moyenne au stade de l'échantillon ou calculer et intégrer les paramètres.

De plus, même après une seule modélisation, les caractéristiques du système peuvent changer (dérive). Le défi est de savoir comment s'adapter à ce moment-là, et si cela peut être réalisé sur un petit programme pour adapter les paramètres.

En dehors de cela, l'approche consistant à considérer l'ensemble comme un système probabiliste peut également être meilleure. Je voudrais créer un nouveau modèle probabiliste et utiliser MCMC etc. pour le trouver comme prochain problème à examiner.

De côté

«Comme j'étais occupé il y a quelques années, j'ai demandé à un collègue de trouver une formule approximative et de la résoudre. ――Mais bien que la formule d'approximation ait été créée, un autre collègue cherchait des paramètres à la main (un dilemme que je ne pouvais pas prendre le temps de résoudre par moi-même ...). J'ai pu l'approcher comme elle était, et elle a été laissée telle quelle parce qu'elle n'était pas particulièrement critique. ―― En fait, la formule d'approximation n'utilisait que le terme différentiel (équivalent au modèle 3 cette fois). Je n'étais pas sceptique, mais j'étais très sceptique quant à savoir si c'était le meilleur. «C'est cette activité que j'ai essayé d'approfondir dans une certaine mesure, y compris la raison pour laquelle je n'ai pas utilisé l'intégration.

À la suite d'une enquête persistante cette fois, il semble que ce qui suit peut être utilisé.

• Il est désormais possible d'automatiser l'optimisation des paramètres.
• J'ai trouvé au moins une approximation qui peut être mise en œuvre dans un système meilleur et plus petit.