### General Shock Model Program

I attached source code for a generalized command-line driven program that takes input discovery data (+ other parameters) and which then generates production curves from the date of first discovery to the present time (and on out into the future). The code basically refactors the specific code from other posts in the Oil Shock Model series, to get a better separation of data from the algorithms.

The program reads the input parameters from either a set of files (for yearly data) or from the term following the option (for fixed constants). It doesn't matter what units you use as long as they remain consistent with the other ones. The program sends the output data to the shell terminal.

- -d
- File containing discovery data (continuous years, one entry per line)
- -s
- File containing shock data (only years with discontinuities), two entries per line, "year" + (space) + "extraction rate that year". The program interpolates between entries, so to prevent shocks, give it only two lines (start and finish times) with the same extraction rate on each line.
- -p
- Optional file containing production data (continuous years, one entry per line). The program doesn't do anything with this data; it only carries it over to make it convenient for comparison to the estimated production.
- -1
- The 1/e rate for remaining fallow (i.e. 1 over the average time between discovery and construction)
- -2
- The 1/e rate for finishing construction (i.e. 1 over the build time)
- -3
- The 1/e rate for reaching maturing (i.e. 1 over the startup time)
- -t
- Timespan for complete simulation

Typical invocation:

`shockmodel.exe -1 0.2 -2 0.2 -3 0.2 -t 600 -d strikes.dat -p prod.dat -s shocks.dat > output`

with Text_IO;

with Ada.Numerics.Elementary_Functions;

with GNAT.Command_Line;

with Ada.Float_Text_IO;

procedure ShockModel is

type Flt is array (Natural range <>) of Float;

type Flt_Access is access all Flt;

Expansion_Factor : constant := 10; -- Year subtics

Time_Span : Integer := 600;

-- Safe array retrieval function

function Get (Arr : Flt; I : Integer) return Float is

begin

if I in Arr'Range then

return Arr (I);

else

return 0.0;

end if;

end Get;

-- Discretized convolution function

function "*" (A, B : in Flt) return Flt is

Total : constant Natural := A'Length + B'Length;

C : Flt (0 .. Total);

V : Float;

begin

for J in 0 .. Total loop

V := 0.0;

for I in 0 .. J loop

V := V + Get (A, I) * Get (B, J - I);

end loop;

C (J) := V;

end loop;

return C;

end "*";

-- Discretized exponential function

function Exponential (L : Natural; Alpha : Float) return Flt is

use Ada.Numerics.Elementary_Functions;

R : Flt (0 .. L - 1);

begin

for I in 0 .. L - 1 loop

R (I) := Alpha * Exp (-Alpha * Float (I));

end loop;

return R;

end Exponential;

-- Fills in data between discrete years

function Expander (Arr : Flt) return Flt is

R : Flt (0 .. Arr'Length * Expansion_Factor - 1);

begin

for I in R'Range loop

R (I) := Arr (I / Expansion_Factor);

end loop;

return R;

end Expander;

-- Safe function for reading forcing function

function Discovery_Window

(Discovery : Flt;

Year : Float;

Start : Float)

return Float

is

Y : constant Float := Year - Start;

begin

if Y < 0.0 then

return 0.0;

end if;

return Discovery (Integer (Y) * Expansion_Factor); -- Subdivide by factor

end Discovery_Window;

type Shock is record

Year : Float; -- Year of shock

Rate : Float; -- Rate of extraction

end record;

type Shocks is array (Natural range <>) of Shock;

type Shocks_Access is access all Shocks;

-- Make the shock a continuous function

-- by interpolating between points in the list

function Interpolate_Shocks (Year : Float; S : Shocks) return Float is

K : Integer := S'Last - 1;

begin

for J in S'First + 1 .. S'Last loop

if S (J).Year > Year then

K := J - 1;

exit;

end if;

end loop;

return S (K).Rate +

(S (K + 1).Rate - S (K).Rate) * (Year - S (K).Year) /

(S (K + 1).Year - S (K).Year);

end Interpolate_Shocks;

-- Strikes into Mature Discovery

Strikes : Flt_Access;

Discovery : Flt_Access;

Phases : Flt (1 .. 3) := (others => 0.0);

-- Generates the mature discovery by repeated convolution

procedure Phasing is

Avg_Fallow : constant Flt :=

Exponential (Time_Span, Phases (1) / Float (Expansion_Factor));

Avg_Build : constant Flt :=

Exponential (Time_Span, Phases (2) / Float (Expansion_Factor));

Avg_Mature : constant Flt :=

Exponential (Time_Span, Phases (3) / Float (Expansion_Factor));

begin

Discovery :=

new Flt'(Expander (Strikes.all) * Avg_Fallow * Avg_Build * Avg_Mature);

end Phasing;

S : Shocks_Access;

Production : Flt_Access; -- To compare against, optional

procedure Compute is

C : Float := 0.0; -- Reserve

Rate : Float; -- Extraction Rate

P : Float; -- Production Rate

DT : constant Float := 0.01; -- 1/100th of year

Start : constant Float := S (S'First).Year;

Finish : constant Float := S (S'Last).Year;

Time : Float := Start;

V : Float := 0.0; -- Volume of oil extracted

Data : Float := 0.0;

PIndex : Integer;

begin

loop

-- Integration of discovery window against extraction rate

Rate := Interpolate_Shocks (Time, S.all);

C := C +

Discovery_Window (Discovery.all, Time, Start) * DT -

C * Rate * DT;

P := Rate * C;

V := V + Rate * C * DT;

if Time > Float (Integer (Time)) - DT / 2.0 and

Time < Float (Integer (Time)) + DT / 2.0

then -- Print output only once per year

PIndex := Integer (Time - Start);

if PIndex in Production'Range then

Data := Production (PIndex);

else

Data := 0.0;

end if;

Text_IO.Put_Line

(Integer'Image (Integer (Time)) & "," & P'Img & "," & Data'Img);

end if;

Time := Time + DT;

exit when Time > Finish;

end loop;

end Compute;

function Get_Flt_Data (File : in String) return Flt_Access is

FT : Text_IO.File_Type;

function Get_Data (Data : in Flt) return Flt is

Value : Flt (0 .. 0);

begin

Ada.Float_Text_IO.Get (FT, Value (0));

Text_IO.Skip_Line (FT);

return Get_Data (Data & Value);

exception

when Text_IO.End_Error =>

return Data;

end Get_Data;

begin

Text_IO.Open (FT, Text_IO.In_File, File);

return new Flt'(Get_Data (Flt'(1 .. 0 => 0.0)));

end Get_Flt_Data;

function Get_Shock_Data (File : in String) return Shocks_Access is

FT : Text_IO.File_Type;

function Get_Data (Data : in Shocks) return Shocks is

Pair : Shocks (0 .. 0);

begin

Ada.Float_Text_IO.Get (FT, Pair (0).Year);

Ada.Float_Text_IO.Get (FT, Pair (0).Rate);

Text_IO.Skip_Line (FT);

return Get_Data (Data & Pair);

exception

when Text_IO.End_Error =>

return Data;

end Get_Data;

begin

Text_IO.Open (FT, Text_IO.In_File, File);

return new Shocks'(Get_Data (Shocks'(1 .. 0 => (0.0, 0.0))));

end Get_Shock_Data;

begin

-- Command line options

loop

case GNAT.Command_Line.Getopt ("1: 2: 3: t: d: s: p:") is

when ASCII.NUL =>

exit;

when '1' =>

Phases (1) := Float'Value (GNAT.Command_Line.Parameter);

when '2' =>

Phases (2) := Float'Value (GNAT.Command_Line.Parameter);

when '3' =>

Phases (3) := Float'Value (GNAT.Command_Line.Parameter);

when 't' =>

Time_Span := Integer'Value (GNAT.Command_Line.Parameter);

when 'd' =>

Strikes := Get_Flt_Data (GNAT.Command_Line.Parameter);

when 's' =>

S := Get_Shock_Data (GNAT.Command_Line.Parameter);

when 'p' =>

Production := Get_Flt_Data (GNAT.Command_Line.Parameter);

when others =>

raise Program_Error; -- cannot occur!

end case;

end loop;

Phasing;

Compute;

exception

when GNAT.Command_Line.Invalid_Switch =>

Text_IO.Put_Line ("Invalid Switch " & GNAT.Command_Line.Full_Switch);

when GNAT.Command_Line.Invalid_Parameter =>

Text_IO.Put_Line ("No parameter for " & GNAT.Command_Line.Full_Switch);

end ShockModel;

## 4 Comments:

Sweet! thank you!

You are most welcome.

I've implemented your code in Matlab and R. I've tried to validate my implementation using your example on ASPO data:

http://mobjectivist.blogspot.com/2005/10/shock-model-with-aspo-discovery-data.html

I retrieved your result only if:

Expansion_Factor= 1;

if Expansion_Factor= 10, the shock model is giving smaller values by 10% compared to the real data.

Khebab,

The "expansion factor" is only used to get better time resolution for high rates. It's a kind of Nyquist factor hack for subdividing a data set into subyear increments.

The best way to validate the porting of the code is to take a delta discovery value (K*delta(t)) and make all the rates the same (the 3 stages plus the depletion rate). The resultant curve should lay exactly on top of a Gamma distribution of order N=4. This is essentially an exponential convolved 3 times with itself. This should show an analytical equivalence.

I mention this because it well tell you if something got screwed up in the expansion factor algorithm code translation. Otherwise, keeping the factor at 1 is OK as long as the numeric errors don't creep in as you increase the rates.

BTW, I also have Python code for this algorithm.

Post a Comment

<< Home

"Like strange bulldogs sniffing each other's butts, you could sense wariness from both sides"