greentangerine
New member
- Local time
- Yesterday, 20:48
- Joined
- Oct 18, 2017
- Messages
- 1
Hi All,
Have programmed this code for a water security issue.
The prob.mod was supplied to use along with our own program.
nrepmax needs to run for 500 000 replicates. It works for around 100 000 replicates, but after that there is an error from the prod.mod code. I have highlighted the error message that comes up, not sure if the highlight works.
Thanks in Advance!
Imports System.Math
Module warragamba
End Module
Have programmed this code for a water security issue.
The prob.mod was supplied to use along with our own program.
nrepmax needs to run for 500 000 replicates. It works for around 100 000 replicates, but after that there is an error from the prod.mod code. I have highlighted the error message that comes up, not sure if the highlight works.
Thanks in Advance!
Imports System.Math
Module warragamba
Code:
Sub main()
' declare variables
Dim nrep, nyear, nrepmax As Integer ' loop variables
Dim demand, volume, capacity, area, evap, VolumeNew As Double
Dim restrict, desal, empty As Integer ' counters
Dim Zt, Z_old, lambda, mean, stdev, correction As Double
Dim percentfull As Double
Dim Prob_restrict, Prob_desal, Prob_empty As Double
' initial values
demand = 400000
capacity = 1886000
nrepmax = [COLOR="red"]100[/COLOR]
'random variable values
lambda = 0.1
mean = 29.277
stdev = 3.589
correction = 0.2184
' replicate loop
For nrep = 1 To nrepmax
' initial loop values
volume = 0.5 * capacity
restrict = 0
desal = 0
empty = 0
' year loop
For nyear = 1 To 40
' increase demand each year (population increase)
demand = demand ^ (1.02)
' restrictions check
If volume < 0.45 * capacity Then
demand = 0.75 * 400000
restrict = restrict + 1
Else
demand = 400000
End If
' Find surface area
area = 17.1 + 57.9 * (volume / 1886000) ^ 0.7555 'm2
' Find evaporation volume
evap = (1530 / 1000) * 0.7 * area ' m3
' Find annual stream inflow into warragamba
Zt = mean + correction * (Z_old - mean) + normal(0.0, stdev)
Qinflow = (Zt * lambda + 1) ^ (1 / lambda)
Zold = Zt ' for next loop
' Calculate volume of reservoir
VolumeNew = volume - demand - evap + Qinflow
' Desalination plant check
If VolumeNew <= 0.3 * capacity Then
'desalination plant built in 2 years
For nyear = nyear + 2 To 40
VolumeNew = VolumeNew + (250 * 365)
Next nyear
desal = desal + 1
End If
' Make sure vol < cap
If VolumeNew > capacity Then
VolumeNew = capacity
ElseIf VolumeNew <= capacity Then
VolumeNew = VolumeNew
End If
' Find % full
percentfull = (VolumeNew / capacity) * 100
' count # times empty
If percentfull = 0 Then
empty = empty + 1
End If
' update volume for next loop
volume = VolumeNew
Next nyear
Next nrep
' calculate probabilities
' restrictions
Prob_restrict = (restrict / CDbl(nrepmax)) * 100
' desalination
Prob_desal = (desal / CDbl(nrepmax)) * 100
' empty
Prob_empty = (empty / CDbl(nrepmax)) * 100
' output
' probability of empty
If Prob_empty > 0.1 Then
Console.WriteLine("The probability of Warragamba Dam being empty is: " & Prob_empty)
Console.WriteLine("Therefore, it will be necessary to augment the system with a new low cost source of water")
Else
Console.WriteLine("The probability of Warragamba Dam being empty is: " & Prob_empty)
Console.WriteLine("The water supply is sufficient")
End If
'probability of restrictions
Console.WriteLine("The probability of Sydney having water restrictions is: " & Prob_restrict)
'probability of desal
Console.WriteLine("The probability of Sydney needing a Desalination Plant to be built is: " & Prob_restrict)
Console.ReadKey()
End Sub
End Module
Module probMod
'
' Module for sampling from probability distributions
' The following module procedures are public:
' uniRan --> uniRan ()
' Function returns a random number sampled from a uniform distribution with limits 0 to 1
' normal --> normal (mean as double, stdDev as double) as double
' Function returns a random number or deviate sampled from a normal
' distribution with Expected value mean and standard devaition stdDev
' poidev --> poidev (mean as double) as long
' Function returns a random number or deviate sampled from a Possion
' distribution with expected value equal to the variable mean
' createHistTable --> CreateHistTable(x() As Double, nX As Long, nBins As Long, firstBin As Double, _
' incrBin As Double, densHist As Boolean, histBins() As String, _
' binCounts() As Double, errorFlag As Boolean, errorMessage As String)
' Subroutine generates a histogram table which may be printed in the calling
' code. The sub supports frequency histograms and density histograms. See agrument
' list below:
' INPUTS:
' x = samples to construct histogram from
' Double type array ()
' nX = number of samples to use in histogram calcs
' Long type scalar
' nBins = number of histogram bins (columns/bars in plot)
' Long type scalar
' **** NOTE: histogram bin and bin count lists are of length nBins+1,
' i.e. have max index nBins and min index 0
' firstBin = smallest value in historagm bin list
' Double type scalar
' incrBin = bin interval or width
' Double type scalar
' densHist = create density historgram if true
' Boolean type
' OUTPUTS:
' histBins = list of histogram bin labels
' String type array
' binCounts = returns frequency or density of each histogram bin
' Double type array
' errorFlag = Boolean type, TRUE if error has occured
' errorMessage = string giving warning or error message
Function normal(mean As Double, stdDev As Double) As Double
'----
'
If stdDev <= 0.0 Then
Call fatalError("f-normal/Negative standard deviation")
Else
normal = mean + stdDev * ppnd(uniRan())
End If
End Function
Function ppnd(p As Double) As Double
'-----
' Algorithm AS 111 Applied Statistics (1977) vol 26, no.1
' "The percentage points of the normal distribution" by J.D. Beasley and S.G. Springer
Dim q, r As Double
Const zero As Double = 0.0, split = 0.42, half = 0.5, one = 1.0,
a0 = 2.50662823884, a1 = -18.61500062529,
a2 = 41.39119773534, a3 = -25.44106049637,
b1 = -8.4735109309, b2 = 23.08336743743,
b3 = -21.06224101826, b4 = 3.13082909833,
c0 = -2.78718931138, c1 = -2.29796479134,
c2 = 4.85014127135, c3 = 2.32121276858,
d1 = 3.54388924762, d2 = 1.63706781897
'-----
q = p - half
If Abs(q) <= split Then
r = q * q
ppnd = q * (((a3 * r + a2) * r + a1) * r + a0) / ((((b4 * r + b3) * r + b2) * r + b1) * r + one)
Else
r = p
If q > zero Then
r = one - p
End If
If r <= zero Then
Call fatalError[COLOR="Red"]("f-ppnd/probability argument p is not in range 0 to 1")[/COLOR]
End If
r = Sqrt(-Log(r))
ppnd = (((c3 * r + c2) * r + c1) * r + c0) / ((d2 * r + d1) * r + one)
If q < zero Then
ppnd = -ppnd
End If
End If
End Function
Function uniRan() As Double
uniRan = Rnd()
End Function
Function poidev(xm As Double) As Long
'------
' Poisson random number generator
'
' Based on code in 'Numerical Recipes' by Press et al. (1986)
'----
Dim pi As Double = 3.141592654
Dim em, t, y As Double
Dim alxm, g, sq As Double
'------
If xm < 12.0 Then
g = Exp(-xm)
em = -1
t = 1.0
Do
em = em + 1.0
t = t * uniRan()
If t <= g Then
Exit Do
End If
Loop
Else
sq = Sqrt(2.0 * xm)
alxm = Log(xm)
g = xm * alxm - gammln(xm + 1.0)
Do
y = Tan(pi * uniRan())
em = sq * y + xm
If em < 0.0 Then
Continue Do
End If
em = Int(em)
t = 0.9 * (1.0 + y ^ 2) * Exp(em * alxm - gammln(em + 1.0) - g)
If uniRan() <= t Then
Exit Do
End If
Loop
End If
poidev = Round(em)
End Function
Sub fatalError(text As String)
Console.WriteLine(text)
Stop
End Sub
Function gammln(xx As Double) As Double
Dim j As Integer
Dim ser, tmp, x, y As Double
Dim cof() As Double = {76.180091729471457, -86.505320329416776,
24.014098240830911, -1.231739572450155,
0.001208650973866179, -0.000005395239384953}
Const stp As Double = 2.5066282746310007
'-----
x = xx
y = x
tmp = x + 5.5
tmp = (x + 0.5) * Log(tmp) - tmp
ser = 1.0000000001900149
For j = 0 To 5
y = y + 1.0
ser = ser + cof(j) / y
Next j
gammln = tmp + Log(stp * ser / x)
End Function
' This subroutine generates a histogram table which may be printed in the calling
' code. The sub supports frequency histograms and density histograms. See agrument
' list below:
' INPUTS:
' x = samples to construct histogram from
' Double type array ()
' nX = number of samples to use in histogram calcs
' Long type scalar
' nBins = number of histogram bins (columns/bars in plot)
' Long type scalar
' **** NOTE: histogram bin and bin count lists are of length nBins+1,
' i.e. have max index nBins and min index 0
' firstBin = smallest value in historagm bin list
' Double type scalar
' incrBin = bin interval or width
' Double type scalar
' densHist = create density historgram if true
' Boolean type
' OUTPUTS:
' histBins = list of histogram bin labels
' String type array
' binCounts = returns frequency or density of each histogram bin
' Double type array
' errorFlag = Boolean type, TRUE if error has occured
' errorMessage = string giving warning or error message
Sub CreateHistTable(x() As Double, nX As Long, nBins As Long, firstBin As Double,
incrBin As Double, densHist As Boolean, histBins() As String,
binCounts() As Double, errorFlag As Boolean, errorMessage As String)
' declare local args
Dim minX As Double
Dim tempBins() As Double
Dim i As Long, j As Long
' initialise
errorFlag = False
' check nBins > 0
If (nBins <= 0) Then
errorFlag = True
errorMessage = "HISTTABLE: ERROR nBins must be > 0"
Exit Sub
End If
' check if minX <= firstBin-incrBin and warn if this is the case
minX = x.Min
If minX <= (firstBin - incrBin) Then
errorMessage = "HISTTABLE: WARNING min(x) is outside lowest bin; histogram may be distorted"
End If
' set up bin values and labels
ReDim tempBins(nBins - 1)
tempBins(0) = firstBin
histBins(0) = CStr(firstBin)
For i = 1 To nBins - 1
tempBins(i) = tempBins(i - 1) + incrBin
histBins(i) = CStr(tempBins(i - 1) + incrBin)
Next i
'histBins(nBins) = "> " & Format(tempBins(nBins - 1), "######0.0#")
histBins(nBins) = "> " & CStr(tempBins(nBins - 1))
For i = 0 To nBins
binCounts(i) = 0
Next i
'--------- histogram calcs
' bin counts
For i = 0 To nX - 1
For j = 0 To nBins - 1
If (x(i) <= tempBins(j)) Then ' update bin count
binCounts(j) = binCounts(j) + 1
Exit For
ElseIf (j = (nBins - 1)) Then ' must belong in > bin
binCounts(nBins) = binCounts(nBins) + 1
End If
Next j
Next i
errorMessage = "HISTTABLE: FrequencyCompletedOK"
' determine density if needed
If densHist Then
For j = 0 To nBins
binCounts(j) = binCounts(j) / (CDbl(nX) * incrBin)
Next j
errorMessage = "HISTTABLE: DensityCompletedOK"
End If
End Sub
Last edited by a moderator:
