Подгонка двухслойной модели для данных профиля ветра в python

Я пытаюсь подогнать модель под свой набор данных профилей ветра, то есть значения скорости ветра u(z) на разных высотах z.

Модель состоит из двух частей, которые я пока упростил до:

u(z) = ust/k * ln(z/z0)     for z < zsl

u(z) = a*z + b              for z > zsl

В логарифмической модели ust и z0 являются свободными параметрами, а k фиксированными. zsl — высота поверхностного слоя, которая также неизвестна априори.

Я хочу подогнать эту модель к своим данным, и я уже пробовал разные подходы. Лучший результат, который я получаю до сих пор, это:

def two_layer(z,hsl,ust,z0,a,b):
    return ust/0.4*(np.log(z/z0)) if z<hsl else a*z+b

two_layer = np.vectorize(two_layer)

def two_layer_err(p,z,u):
    return two_layer(z,*p)-u

popt, pcov ,infodict, mesg, ier = optimize.leastsq(two_layer_err,[150.,0.3,0.002,0.1,10.],(wspd_hgt,prof),full_output=1)

# wspd_hgt are my measurements heights and 
# prof are the corresponding wind speed values

Это дает мне разумные оценки для всех параметров, кроме zsl, который не изменяется во время процедуры подгонки. Я предполагаю, что это связано с тем, что используется как порог, а не как параметр функции. Можно ли как-то изменить zsl во время оптимизации?

Я пробовал что-то с numpy.piece, но это не сработало, возможно, потому, что я не очень хорошо это понимаю, или я мог бы быть совершенно не здесь, потому что это не подходит для моего дела.

По идее, профиль ветра выглядит так, если оси повернуть вспять (на графике z и u): График профиля ветра


python model-fitting
person Peter9192    schedule 04.03.2016    source источник
comment
Я думаю, ваша проблема в том, что функция в общем случае не является непрерывной, поэтому при определении производной по в.р.т. hsl/zsl вы, вероятно, получите бесконечность, а curve_fit не справится с этим. Я замечал, что подобные вещи происходили и раньше при добавлении параметров, поведение которых не является «математическим», но имеет смысл с точки зрения программирования. Я попытался поиграть с вашей проблемой с другим пакетом, но пока не повезло. Есть ли способ заставить функцию быть непрерывной?   -  person tBuLi    schedule 26.03.2016
comment
@tBuLi, спасибо за ваш комментарий. К сожалению, нет способа сделать это. Лучшее, что я придумал до сих пор, - это перебрать разумные значения hsl и подогнать логарифмическую функцию к точкам ниже и линейной функции выше.   -  person Peter9192    schedule 29.03.2016

Ответы (1)


arrow_upward
1
arrow_downward

Я думаю, что у меня наконец-то есть решение для этого типа проблемы, с которым я столкнулся, отвечая на похожий вопрос.

Решение, по-видимому, состоит в том, чтобы реализовать ограничение, говорящее, что u1 == u2 при переключении между двумя моделями. Поскольку я не могу попробовать это с вашей моделью из-за отсутствия опубликованных данных, вместо этого я покажу, как это работает для другой модели, и вы можете адаптировать ее к своей ситуации. Я решил это, используя оболочку scipy, которую я написал, чтобы сделать решение таких проблем более питоническим, под названием symfit. Но вы можете сделать то же самое, используя алгоритм SLSQP в scipy, если хотите.

from symfit import parameters, variables, Fit, Piecewise, exp, Eq
import numpy as np
import matplotlib.pyplot as plt

t, y = variables('t, y')
m, c, d, k, t0 = parameters('m, c, d, k, t0')

# Help the fit by bounding the switchpoint between the models
t0.min = 0.6
t0.max = 0.9

# Make a piecewise model
y1 = m * t + c
y2 = d * exp(- k * t)
model = {y: Piecewise((y1, t <= t0), (y2, t > t0))}

# As a constraint, we demand equality between the two models at the point t0
# Substitutes t in the components by t0
constraints = [Eq(y1.subs({t: t0}), y2.subs({t: t0}))]

# Read the data
tdata, ydata = np.genfromtxt('Experimental Data.csv', delimiter=',', skip_header=1).T

fit = Fit(model, t=tdata, y=ydata, constraints=constraints)
fit_result = fit.execute()
print(fit_result)

plt.scatter(tdata, ydata)
plt.plot(tdata, fit.model(t=tdata, **fit_result.params).y)
plt.show()

введите здесь описание изображения

Я думаю, вы должны быть в состоянии адаптировать этот пример к вашей ситуации!

Редактировать: как указано в комментарии, можно также потребовать сопоставления производных. Чтобы сделать это, для приведенного выше примера необходим следующий дополнительный код:

from symfit import Derivative

dy1dt = Derivative(y1, t)
dy2dt = Derivative(y2, t)
constraints = [
    Eq(y1.subs({t: t0}), y2.subs({t: t0})),
    Eq(dy1dt.subs({t: t0}), dy2dt.subs({t: t0}))
]

Это должно делать свое дело! Так что с точки зрения программирования это очень выполнимо, хотя в зависимости от модели это может не иметь решения.

person tBuLi    schedule 27.08.2018
comment
У меня нет времени вникать в это прямо сейчас, но можно ли добавить второе ограничение, указывающее, что производные тоже должны совпадать? - person Peter9192; 28.08.2018
comment
@ Peter9192 Peter9192 Да, я обновил свой ответ, чтобы включить это. - person tBuLi; 28.08.2018