How to find the parameters of Herschel-Bulkley model through Nonlinear Regression!


(Maamar Dliouah) #1

Dear all,

First, I would like to thank everyone involved in this magnificent project, Math.NET saved my life!
I have few questions about the linear and nonlinear regression, I am a civil engineer and when I was working on my Master’s degree, I needed to develop a C# application that calculates the Rheological parameters of concrete based on data acquired from a test.
One of the models that describes the rheological behavior of concrete is the “Herschel-Bulkley model” and it has this formula :

y = T + K*x^n

x (the shear-rate), y (shear-stress) are the values obtained from the test, while T,K and N are the parameters I need to determine.
I know that the value of “T” is between 0 and Ymin (Ymin is the smallest data point from the test), so here is what I did:
Since it is nonlinear equation, I had to make it linear, like this : ln(y-T) = ln(K) + n*ln(x)
creat an array of possible values of T, from 0 to Ymin, and try each value in the equation, then through linear regression I find the values of K and N, then calculate the SSD, and store the results in an array, after I finish all the possible values of T, I see which one had the smallest SSD, and use it to find the K and N . this method works, but I feel it is not as smart or elegant as it should be, there must be a better way to do it, and I was hoping to find it here, it is also very slow.
here is the code that I used:

public static double HerschelBulkley(double shearRate, double tau0, double k, double n)
{
    var t = tau0 + k * Math.Pow(shearRate, n);
    return t;
}

public static (double Tau0, double K, double N, double DeltaMin, double RSquared) HerschelBulkleyModel(double[] shear, double[] shearRate, double step = 1000.0)
{
    // Calculate the number values from 0.0 to Shear.Min;
    var sm = (int) Math.Floor(shear.Min() * step);
    // Populate the Array of Tau0 with the values from 0 to sm
    var tau0Array = Enumerable.Range(0, sm).Select(t => t / step).ToArray();
    var kArray = new double[sm];
    var nArray = new double[sm];
    var deltaArray = new double[sm];
    var rSquaredArray = new double[sm];
    var shearRateLn = shearRate.Select(s => Math.Log(s)).ToArray();
    for (var i = 0; i < sm; i++)
    {
        var shearLn = shear.Select(s => Math.Log(s - tau0Array[i])).ToArray();
        var param = Fit.Line(shearRateLn, shearLn);
        kArray[i] = Math.Exp(param.Item1);
        nArray[i] = param.Item2;
        var shearHerschel = shearRate.Select(sr => HerschelBulkley(sr, tau0Array[i], kArray[i], nArray[i])).ToArray();
        deltaArray[i] = Distance.SSD(shearHerschel, shear);
        rSquaredArray[i] = GoodnessOfFit.RSquared(shearHerschel, shear);
    }
    var deltaMin = deltaArray.Min();
    var index = Array.IndexOf(deltaArray, deltaMin);
    var tau0 = tau0Array[index];
    var k = kArray[index];
    var n = nArray[index];
    var rSquared = rSquaredArray[index];
    return (tau0, k, n, deltaMin, rSquared);
}

(Peter Vanderwaart) #2

I know this is an older post, but I just came across it recently, and I have a few comments. As it happens, I have posted two programs for fitting an equivalent equation in the “Exponential Fit” thread. Like you, I noticed that if one of the unknowns if fixed, the other two can be determined by linear regression, so the whole problem can be regarded as an optimization problem in one variable. It your case, that variable was T. In my program it was the equivalent of N.

I think your approach is probably better. The most difficult part of the problem is to distinguish between linear growth created by K and exponential grown created by N. Your approach handles this more directly, and you can get t-statistics and other stats from the regression to assess how good the fit might be.

I never considered using the SSD from the linear regression as the objective function. Great idea! Quite likely better than using the SSD from estimating the original equation as I did.

My only real contribution is to note that if you create a function f(K) which returns the SSD of the regression, then you can use a scalar mimimizer to find the optimum. I used GoldenSectionMinimizer.Minimum().