2

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.piecewise 尝试了一些东西，但效果不太好，可能是因为我不太了解它，或者我可能完全离开这里，因为它不适合我的事业。

对于这个想法，如果轴反转（z绘制对比u），风廓线看起来像这样： 风廓线图

4

1 回答 1

1

I think I finally have a solution for this type of problem, which I came across while answering a similar question.

The solution seems to be to implement a constraint saying that u1 == u2 at the switch between the two models. Because I cannot try it with your model because there is no data posted, I will instead show how it works for the other model, and you can adapt it to your situation. I solved this using a scipy wrapper I wrote to make fitting of such problems more pythonic, called symfit. But you could do the same using the SLSQP algorithm in scipy if you prefer.

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()

enter image description here

I think you should be able to adapt this example to your situation!

Edit: as requested in the comment, it is possible to demand matching derivatives as well. In order to do this, the following additional code is needed for the above example:

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}))
]

That should do the trick! So from a programming point of view it is very doable, although depending on the model this might not actually have a solution.

于 2018-08-27T09:22:22.857 回答