From 9f0b9d28405b7900aa7c6e8a27658c97d80360d7 Mon Sep 17 00:00:00 2001 From: Magne Sjaastad Date: Tue, 13 Oct 2026 10:14:07 +0200 Subject: [PATCH] #11172: Move decline curve math to RigDeclineCurveCalculator Add least squares regression as an alternative to the start/end point fit. --- .../RigDeclineCurveCalculator.cpp | 316 ++++++++++++++++++ .../RigDeclineCurveCalculator.h | 69 ++++ .../RigDeclineCurveCalculator-Test.cpp | 69 ++++ 3 files changed, 454 insertions(+) diff --git a/ApplicationLibCode/ReservoirDataModel/RigDeclineCurveCalculator.cpp b/ApplicationLibCode/ReservoirDataModel/RigDeclineCurveCalculator.cpp index 2bbd7d789e..13dcec6a23 100644 --- a/ApplicationLibCode/ReservoirDataModel/RigDeclineCurveCalculator.cpp +++ b/ApplicationLibCode/ReservoirDataModel/RigDeclineCurveCalculator.cpp @@ -18,7 +18,12 @@ #include "RigDeclineCurveCalculator.h" +#include +#include + +#include #include +#include //-------------------------------------------------------------------------------------------------- /// @@ -95,3 +100,314 @@ double RigDeclineCurveCalculator::computeCumulativeProductionHyperbolicDecline( ( std::pow( initialProductionRateQi, 1.0 - hyperbolicDeclineConstantB ) - std::pow( productionRate, 1.0 - hyperbolicDeclineConstantB ) ); } + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::optional> RigDeclineCurveCalculator::fitLine( const std::vector& t, const std::vector& y ) +{ + if ( t.size() != y.size() || t.size() < 2 ) return std::nullopt; + + // Scale time to avoid numerical problems when time is given in seconds + const double timeScale = 365.25 * 24.0 * 3600.0; + + const double n = static_cast( t.size() ); + double sumT = 0.0, sumY = 0.0, sumTT = 0.0, sumTY = 0.0; + for ( size_t i = 0; i < t.size(); i++ ) + { + const double ts = t[i] / timeScale; + sumT += ts; + sumY += y[i]; + sumTT += ts * ts; + sumTY += ts * y[i]; + } + + const double denominator = n * sumTT - sumT * sumT; + if ( std::abs( denominator ) < 1.0e-12 ) return std::nullopt; + + const double slopePerYear = ( n * sumTY - sumT * sumY ) / denominator; + const double intercept = ( sumY - slopePerYear * sumT ) / n; + if ( !std::isfinite( slopePerYear ) || !std::isfinite( intercept ) ) return std::nullopt; + + return std::make_pair( intercept, slopePerYear / timeScale ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +RigDeclineCurveCalculator::FitResult RigDeclineCurveCalculator::fitExponentialDecline( const std::vector& times, + const std::vector& flowRates ) +{ + // ln(q) = ln(qi) - Di * t + std::vector t, y; + for ( size_t i = 0; i < std::min( times.size(), flowRates.size() ); i++ ) + { + if ( flowRates[i] > 0.0 ) + { + t.push_back( times[i] ); + y.push_back( std::log( flowRates[i] ) ); + } + } + + auto line = fitLine( t, y ); + if ( !line ) return std::nullopt; + + return std::make_pair( std::exp( line->first ), -line->second ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +RigDeclineCurveCalculator::FitResult RigDeclineCurveCalculator::fitHarmonicDecline( const std::vector& times, + const std::vector& flowRates ) +{ + // 1/q = 1/qi + (Di/qi) * t + std::vector t, y; + for ( size_t i = 0; i < std::min( times.size(), flowRates.size() ); i++ ) + { + if ( flowRates[i] > 0.0 ) + { + t.push_back( times[i] ); + y.push_back( 1.0 / flowRates[i] ); + } + } + + auto line = fitLine( t, y ); + if ( !line || line->first <= 0.0 ) return std::nullopt; + + const double qi = 1.0 / line->first; + return std::make_pair( qi, line->second * qi ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +RigDeclineCurveCalculator::FitResult RigDeclineCurveCalculator::fitHyperbolicDecline( const std::vector& times, + const std::vector& flowRates, + double hyperbolicDeclineConstantB ) +{ + if ( hyperbolicDeclineConstantB <= 0.0 ) return std::nullopt; + + // q^-b = qi^-b * (1 + b * Di * t) + std::vector t, y; + for ( size_t i = 0; i < std::min( times.size(), flowRates.size() ); i++ ) + { + if ( flowRates[i] > 0.0 ) + { + t.push_back( times[i] ); + y.push_back( std::pow( flowRates[i], -hyperbolicDeclineConstantB ) ); + } + } + + auto line = fitLine( t, y ); + if ( !line || line->first <= 0.0 ) return std::nullopt; + + const double qi = std::pow( line->first, -1.0 / hyperbolicDeclineConstantB ); + const double di = line->second / ( line->first * hyperbolicDeclineConstantB ); + return std::make_pair( qi, di ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::pair RigDeclineCurveCalculator::computeStartEndPointParameters( const std::vector& values, + const std::vector& timeSteps, + bool isAccumulated ) +{ + auto computeProductionRate = []( double t0, double v0, double t1, double v1 ) { return ( v1 - v0 ) / ( t1 - t0 ); }; + + const double endTime = static_cast( timeSteps.back() ); + + if ( !isAccumulated ) + { + // Use the first (time step, value) pair as t0. The last point (within the selected range) is the initial production rate. + const double initialProductionRate = values.back(); + const double initialDeclineRate = + computeDeclineRate( static_cast( timeSteps.front() ), values.front(), endTime, initialProductionRate ); + return { initialProductionRate, initialDeclineRate }; + } + + // Select a point (t0) 1/4 into the user-specified range + const double historyStep = 0.25; + const size_t idx0 = static_cast( timeSteps.size() * historyStep ); + const double t0 = static_cast( timeSteps[idx0] ); + const double v0 = values[idx0]; + + // For accumulated result: compute the initial production rate from the two points. + const double initialProductionRate = computeProductionRate( t0, v0, endTime, values.back() ); + + // Compute the production rate at time t0 by using the first point in the existing curve. + const double productionRate0 = computeProductionRate( static_cast( timeSteps.front() ), values.front(), t0, v0 ); + + const double initialDeclineRate = computeDeclineRate( t0, productionRate0, endTime, initialProductionRate ); + return { initialProductionRate, initialDeclineRate }; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::optional + RigDeclineCurveCalculator::computeDeclineParameters( const DeclineModel& model, + const std::vector& values, + const std::vector& timeSteps ) +{ + if ( values.empty() || values.size() != timeSteps.size() ) return std::nullopt; + + const time_t referenceTime = timeSteps.back(); + + double initialProductionRate = 0.0; + double initialDeclineRate = 0.0; + + if ( model.fitMethod == FitMethod::START_END_POINT ) + { + std::tie( initialProductionRate, initialDeclineRate ) = computeStartEndPointParameters( values, timeSteps, model.isAccumulated ); + } + else + { + // Regression is performed on flow rates. Derive rates from the differences if the curve is accumulated. + std::vector times; + std::vector rates; + for ( size_t i = 0; i < values.size(); i++ ) + { + if ( model.isAccumulated ) + { + if ( i == 0 || timeSteps[i] <= timeSteps[i - 1] ) continue; + rates.push_back( ( values[i] - values[i - 1] ) / static_cast( timeSteps[i] - timeSteps[i - 1] ) ); + + // The average rate over the interval is best represented at the interval midpoint + times.push_back( 0.5 * static_cast( timeSteps[i - 1] - referenceTime ) + + 0.5 * static_cast( timeSteps[i] - referenceTime ) ); + } + else + { + rates.push_back( values[i] ); + times.push_back( static_cast( timeSteps[i] - referenceTime ) ); + } + } + + FitResult fit; + if ( model.type == DeclineCurveType::EXPONENTIAL ) + { + fit = fitExponentialDecline( times, rates ); + } + else if ( model.type == DeclineCurveType::HARMONIC ) + { + fit = fitHarmonicDecline( times, rates ); + } + else + { + fit = fitHyperbolicDecline( times, rates, model.hyperbolicDeclineConstantB ); + } + + if ( !fit ) return std::nullopt; + std::tie( initialProductionRate, initialDeclineRate ) = *fit; + } + + if ( !std::isfinite( initialProductionRate ) || !std::isfinite( initialDeclineRate ) ) return std::nullopt; + + return DeclineParameters{ initialProductionRate, initialDeclineRate, referenceTime }; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +double RigDeclineCurveCalculator::computePredictedValue( const DeclineModel& model, const DeclineParameters& parameters, double timeSinceStart ) +{ + const double qi = parameters.initialProductionRate; + const double di = parameters.initialDeclineRate; + const double b = model.hyperbolicDeclineConstantB; + + switch ( model.type ) + { + case DeclineCurveType::EXPONENTIAL: + return model.isAccumulated ? computeCumulativeProductionExponentialDecline( qi, di, timeSinceStart ) + : computeFlowRateExponentialDecline( qi, di, timeSinceStart ); + case DeclineCurveType::HARMONIC: + return model.isAccumulated ? computeCumulativeProductionHarmonicDecline( qi, di, timeSinceStart ) + : computeFlowRateHarmonicDecline( qi, di, timeSinceStart ); + case DeclineCurveType::HYPERBOLIC: + return model.isAccumulated ? computeCumulativeProductionHyperbolicDecline( qi, di, timeSinceStart, b ) + : computeFlowRateHyperbolicDecline( qi, di, timeSinceStart, b ); + } + + return 0.0; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::optional> RigDeclineCurveCalculator::computeDeclineCurveValues( const DeclineModel& model, + const std::vector& values, + const std::vector& timeSteps, + const std::vector& futureTimeSteps ) +{ + auto parameters = computeDeclineParameters( model, values, timeSteps ); + if ( !parameters ) return std::nullopt; + + auto computeValue = [&]( time_t time ) + { + double timeSinceStart = static_cast( time - parameters->referenceTime ); + double predictedValue = computePredictedValue( model, *parameters, timeSinceStart ); + if ( model.isAccumulated ) predictedValue += values.back(); + return predictedValue; + }; + + // The history part is a backward extrapolation of the decline curve. Keep the original value if the equation is not valid in the + // history time range. + std::vector outValues; + for ( size_t i = 0; i < values.size(); i++ ) + { + double value = computeValue( timeSteps[i] ); + outValues.push_back( std::isfinite( value ) ? value : values[i] ); + } + + for ( time_t futureTime : futureTimeSteps ) + { + outValues.push_back( computeValue( futureTime ) ); + } + + return outValues; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +QString RigDeclineCurveCalculator::equationText( const DeclineModel& model, const DeclineParameters& parameters ) +{ + QString rateEquation; + QString cumulativeEquation; + switch ( model.type ) + { + case DeclineCurveType::EXPONENTIAL: + rateEquation = "q(t) = qi * exp(-Di * t)"; + cumulativeEquation = "Np(t) = Np0 + (qi - q(t)) / Di"; + break; + case DeclineCurveType::HARMONIC: + rateEquation = "q(t) = qi / (1 + Di * t)"; + cumulativeEquation = "Np(t) = Np0 + (qi / Di) * ln(qi / q(t))"; + break; + case DeclineCurveType::HYPERBOLIC: + rateEquation = "q(t) = qi / (1 + b * Di * t)^(1 / b)"; + cumulativeEquation = "Np(t) = Np0 + qi^b / ((1 - b) * Di) * (qi^(1 - b) - q(t)^(1 - b))"; + break; + } + + QStringList lines; + lines << rateEquation; + if ( model.isAccumulated ) lines << cumulativeEquation; + + // qi is derived per second for accumulated curves, scale to the time unit of t + const double qi = model.isAccumulated ? parameters.initialProductionRate * secondsPerYear : parameters.initialProductionRate; + const QString qiUnit = model.isAccumulated ? " [per year]" : ""; + QString parameterText = + QString( "qi = %1%2, Di = %3 [1/year]" ).arg( qi, 0, 'g', 6 ).arg( qiUnit ).arg( parameters.initialDeclineRate * secondsPerYear, 0, 'g', 6 ); + if ( model.type == DeclineCurveType::HYPERBOLIC ) + parameterText += QString( ", b = %1" ).arg( model.hyperbolicDeclineConstantB, 0, 'g', 4 ); + lines << parameterText; + + const QString referenceDate = QDateTime::fromSecsSinceEpoch( parameters.referenceTime ).date().toString( Qt::ISODate ); + lines << QString( "t is time in years since %1" ).arg( referenceDate ); + if ( model.isAccumulated ) lines << "Np0 is the accumulated value at t = 0"; + + return lines.join( "\n" ); +} diff --git a/ApplicationLibCode/ReservoirDataModel/RigDeclineCurveCalculator.h b/ApplicationLibCode/ReservoirDataModel/RigDeclineCurveCalculator.h index 4a035fbf55..30184024db 100644 --- a/ApplicationLibCode/ReservoirDataModel/RigDeclineCurveCalculator.h +++ b/ApplicationLibCode/ReservoirDataModel/RigDeclineCurveCalculator.h @@ -18,6 +18,13 @@ #pragma once +#include + +#include +#include +#include +#include + //================================================================================================== // // @@ -26,6 +33,51 @@ class RigDeclineCurveCalculator { public: + enum class DeclineCurveType + { + EXPONENTIAL, + HARMONIC, + HYPERBOLIC + }; + + enum class FitMethod + { + START_END_POINT, + REGRESSION + }; + + struct DeclineModel + { + DeclineCurveType type; + FitMethod fitMethod; + double hyperbolicDeclineConstantB; + bool isAccumulated; + }; + + struct DeclineParameters + { + double initialProductionRate; // [rate unit] + double initialDeclineRate; // [1/second] + time_t referenceTime; // time at which the initial production rate applies + }; + + static constexpr double secondsPerYear = 365.25 * 24.0 * 3600.0; + + // Estimate the decline parameters from the history values. The reference time is the last time step. + static std::optional + computeDeclineParameters( const DeclineModel& model, const std::vector& values, const std::vector& timeSteps ); + + // Flow rate, or cumulative production increment since the reference time if the model is accumulated + static double computePredictedValue( const DeclineModel& model, const DeclineParameters& parameters, double timeSinceStart ); + + // Returns the decline curve evaluated at the history time steps (backward extrapolation) followed by the future time steps. + static std::optional> computeDeclineCurveValues( const DeclineModel& model, + const std::vector& values, + const std::vector& timeSteps, + const std::vector& futureTimeSteps ); + + static QString equationText( const DeclineModel& model, const DeclineParameters& parameters ); + static double computeDeclineRate( double time0, double value0, double time1, double value1 ); static double computeFlowRateExponentialDecline( double initialProductionRateQi, double initialDeclineRateDi, double timeSinceStart ); @@ -47,4 +99,21 @@ public: double initialDeclineRateDi, double timeSinceStart, double hyperbolicDeclineConstantB ); + + // Least squares fit of the decline equation to the given flow rates. Time values are seconds relative to the reference time (the + // initial production rate is the rate at time zero). Returns {initialProductionRateQi, initialDeclineRateDi}, or nullopt if the + // data does not allow a fit. Non-positive rates are ignored. + using FitResult = std::optional>; + + static FitResult fitExponentialDecline( const std::vector& times, const std::vector& flowRates ); + static FitResult fitHarmonicDecline( const std::vector& times, const std::vector& flowRates ); + static FitResult + fitHyperbolicDecline( const std::vector& times, const std::vector& flowRates, double hyperbolicDeclineConstantB ); + +private: + static std::pair + computeStartEndPointParameters( const std::vector& values, const std::vector& timeSteps, bool isAccumulated ); + + // Least squares fit of y = a + b * t. Returns {a, b}. + static std::optional> fitLine( const std::vector& t, const std::vector& y ); }; diff --git a/ApplicationLibCode/UnitTests/RigDeclineCurveCalculator-Test.cpp b/ApplicationLibCode/UnitTests/RigDeclineCurveCalculator-Test.cpp index d39b2531d7..3a1e17586f 100644 --- a/ApplicationLibCode/UnitTests/RigDeclineCurveCalculator-Test.cpp +++ b/ApplicationLibCode/UnitTests/RigDeclineCurveCalculator-Test.cpp @@ -20,6 +20,8 @@ #include "RigDeclineCurveCalculator.h" +#include + // Test case from https://web.archive.org/web/20120710104333/http://infohost.nmt.edu/~petro/faculty/Kelly/Deline.pdf // Example Problem 8-1: // "Given that a well has declined from 100 stb/day to 96 stb/day during a one-month period, @@ -55,3 +57,70 @@ TEST( RigDeclineCurveCalculatorTests, exponentialFlowRateExample ) EXPECT_NEAR( q, 61.27, 0.005 ); } + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +TEST( RigDeclineCurveCalculatorTests, fitRecoversParametersFromSyntheticData ) +{ + const double qi = 100.0; + const double di = 0.4 / ( 365.25 * 24.0 * 3600.0 ); // per second + const double b = 0.5; + + std::vector times; + std::vector exponential, harmonic, hyperbolic; + for ( int i = -20; i <= 0; i++ ) + { + const double t = i * 30.0 * 24.0 * 3600.0; + times.push_back( t ); + exponential.push_back( RigDeclineCurveCalculator::computeFlowRateExponentialDecline( qi, di, t ) ); + harmonic.push_back( RigDeclineCurveCalculator::computeFlowRateHarmonicDecline( qi, di, t ) ); + hyperbolic.push_back( RigDeclineCurveCalculator::computeFlowRateHyperbolicDecline( qi, di, t, b ) ); + } + + auto expFit = RigDeclineCurveCalculator::fitExponentialDecline( times, exponential ); + ASSERT_TRUE( expFit.has_value() ); + EXPECT_NEAR( qi, expFit->first, 1.0e-6 ); + EXPECT_NEAR( di, expFit->second, 1.0e-12 ); + + auto harmonicFit = RigDeclineCurveCalculator::fitHarmonicDecline( times, harmonic ); + ASSERT_TRUE( harmonicFit.has_value() ); + EXPECT_NEAR( qi, harmonicFit->first, 1.0e-6 ); + EXPECT_NEAR( di, harmonicFit->second, 1.0e-12 ); + + auto hyperbolicFit = RigDeclineCurveCalculator::fitHyperbolicDecline( times, hyperbolic, b ); + ASSERT_TRUE( hyperbolicFit.has_value() ); + EXPECT_NEAR( qi, hyperbolicFit->first, 1.0e-6 ); + EXPECT_NEAR( di, hyperbolicFit->second, 1.0e-12 ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +TEST( RigDeclineCurveCalculatorTests, fitFailsWithTooFewValidPoints ) +{ + EXPECT_FALSE( RigDeclineCurveCalculator::fitExponentialDecline( { -10.0, 0.0 }, { 0.0, 5.0 } ).has_value() ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +TEST( RigDeclineCurveCalculatorTests, historyIsBackwardExtrapolationOfFittedCurve ) +{ + const double day = 24.0 * 3600.0; + + using Calc = RigDeclineCurveCalculator; + Calc::DeclineModel model{ Calc::DeclineCurveType::EXPONENTIAL, Calc::FitMethod::REGRESSION, 0.5, false }; + + // Noisy history around an exponential decline + std::vector timeSteps = { 0, static_cast( 100 * day ), static_cast( 200 * day ), static_cast( 300 * day ) }; + std::vector values = { 100.0, 80.0, 66.0, 50.0 }; + + auto result = Calc::computeDeclineCurveValues( model, values, timeSteps, { static_cast( 400 * day ) } ); + ASSERT_TRUE( result.has_value() ); + ASSERT_EQ( 5u, result->size() ); + + EXPECT_NE( values[0], ( *result )[0] ); + EXPECT_GT( ( *result )[0], ( *result )[1] ); + EXPECT_GT( ( *result )[3], ( *result )[4] ); +} \ No newline at end of file