拟合两层模型以在python中缠绕轮廓数据 [英] Fitting a two-layer model to wind profile data in python

问题描述

我正在尝试将模型拟合到我的风廓线数据集中，即在不同高度z处的风速值u(z).

I'm trying to fit a model to my dataset of wind profiles, i.e. wind speed values u(z) at different altitudes z.

该模型由两部分组成，我现在将其简化为:

The model consists of two parts, which I for now simplified to:

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

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

在对数模型中，ust和z0是自由参数，k是固定的. zsl是表面层的高度，这也不是先验知识.

In the logarithmic model, ust and z0 are free parameters k is fixed. zsl is the height of the surface layer, which is also not known a priori.

我想将此模型适合我的数据，并且我已经尝试了不同的方法.到目前为止，我得到的最好结果是:

I want to fit this model to my data, and I have already tried different approaches. The best result I'm getting so far is with:

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在优化过程中变化?

This gives me reasonable estimates for all parameters, except for zsl which is not changed during the fitting procedure. I guess this has to do with the fact that is used as a threshold rather than a function parameter. Is there any way I could let zsl be varied during the optimization?

我用numpy.piecewise尝试了一些操作，但是效果不佳，可能是因为我不太了解它，或者我可能完全不在这里，因为它不适合我的事业.

I tried something with numpy.piecewise, but that didn't work very well, perhaps because I don't really understand it very well, or I might be completely off here because it's not suitable for my cause.

对于这个想法，如果轴反转，则风廓形如下(z与u的关系图):

For the idea, the wind profile looks like this if the axes are reversed (z plotted versus u):

推荐答案

我认为我终于有了解决此类问题的解决方案，我在回答

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

解决方案似乎是在两个模型之间的切换处实施约束u1 == u2.因为没有发布数据，所以我无法在您的模型上进行尝试，所以我将展示它在其他模型上的工作方式，您可以根据自己的情况进行调整.我使用编写的scipy包装器解决了这个问题，以使这种问题更符合Pythonic，称为symfit.但是，如果愿意，可以使用scipy中的SLSQP算法进行相同的操作.

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

我认为您应该能够使此示例适应您的情况！

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

根据注释中的要求，也可以要求匹配的派生词.为此，上面的示例需要以下附加代码:

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.

这篇关于拟合两层模型以在python中缠绕轮廓数据的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！