RUMUS YANG DIGUNAKAN
Prediksi pasang surut (pasut) menggunakan rumus Harmonic Analysis dengan metode least square untuk mencari 9 constituents utama pembangkit pasut dan 9 phase.
Dalam referensi tertulis,
"The ideal tide curve for any given port is represented as an average height Z0 plus a sum of terms ("constituents'') each of which is of the form f(t) = H cos(at + \phi). The time t is measured in hours, and f comes out in feet. The numbers H,a, \phi are the amplitude, the speed and the phase of the constituent."
Karena data hasil bacaan pasut yang tersedia adalah dalam meter dan dicatat tiap jam, maka nilai f tidak lagi dinyatakan dalam feet, melainkan dalam meter.
Rumus kemudian dikembangkan menjadi:
dalam format matrix dapat ditulis:
Jika MatrixX telah didapat, maka dengan menggunakan rumus(4) dan rumus(5),amplitudo dan phase tiap constituent dapat dihitung.
INPUT (DATA HASIL PENGUKURAN MUKA AIR)
Hasil pengukuran muka air didapat dari pengamatan selama 29 hari dengan pencatatan tiap jam, dimulai jam 0:00 sampai jam 23:00 tanggal 1 Maret 2008 sampai jam 23:00 tanggal 29 Maret 2008 yang ditabelkan sebagai berikut:
PROSES PEMBACAAN DATA DAN PERHITUNGAN
Proses pembacaan dan perhitungan dengan menggunakan program Visual Basic Application (macro) yang ada di microsoft excel 2003.
Khusus untuk prosedure atau program perkalian matirx dan inverse didapat dari http://www.alglib.net/
Pada file excel spread sheet saya, data pengukuran muka air pertama (1.850m) ada di cell $C$10 sedangkan data terakhir (2.110m) ada di $Z$38, sehingga range data pengukuran adalah di $C$10:$Z$38.
Pastikan bahwa di bagian Declaration di awal program diset :
Option Explicit 'berguna untuk mendeteksi definisi variable
Option Base 1 'hitungan matrix dimulai dari 1, kalo tidak diset, maka default index matrix adalah 0
Sub PrediksiPasutDenganLeastSquare()
'1. membaca data pengukuran muka air
dim cr as Range, rgData as Range
set rgData=Activesheet.Range("$C$10:$Z$38")
For Each cr In rgData.Columns(1).Cells
i = i + 1
j = 1 + (i - 1) * 24
ReDim Preserve MatrixL(1 To j + 23)
For jm = 0 To 23
MatrixL(j + jm) = cr.Offset(0, jm).Value
Next jm
Next
'2. membaca atau set variable periode tiap Constituents
Dim w(1 To 9) As Double 'periode dari 9 Constituents pasut
Dim pi As Double
pi = 4 * Atn(1)
w(1) = 2# * pi / 12.4206 'M2
w(2) = 2# * pi / 12# 'S2
w(3) = 2# * pi / 12.6582 'N2
w(4) = 2# * pi / 11.9673 'K2
w(5) = 2# * pi / 23.9346 'K1
w(6) = 2# * pi / 25.8194 'O1
w(7) = 2# * pi / 24.0658 'P1
w(8) = 2# * pi / 6.2103 'M4
w(9) = 2# * pi / 6.1033 'MS4
'3. membuat atau membetuk MatrixA() atau Matrix coeffisien
Dim it As Integer,MatrixA() as Double
ReDim Preserve matrixA(1 To UBound(MatrixL), 1 To 19)
For i = 1 To UBound(MatrixL)
matrixA(i, 1) = 1
it = i
For j = 1 To 9
matrixA(i, 2 * j) = Cos(w(j) * it)
matrixA(i, 2 * j + 1) = -Sin(w(j) * it)
Next j
Next i
'4. proses least square
Dim weight() As Double: ReDim Preserve weight(UBound(MatrixL, 1))
For i = LBound(weight) To UBound(weight): weight(i) = 1: Next i 'matrix bobot=matrix identitas
MatrixX() = clsLSQ.LSPAR(matrixA, MatrixL, weight)
'clsLSQ adalah class yang saya buat untuk proses perhitungan least square, tidak dibahas di session ini
'5. menghitung amplitudo dan phase 9 konstituent dan sekaligus menampilkan hasilnya di excel
Dim A As Double, B As Double, H(1 To 9) As Double, Phase(1 To 9) As Double
Dim ph As Double, Zo As Double
Const addPrint As String = "H66"
With Range(addPrint)
.Offset(0, 3) = MatrixX(1)'mencetak Zo atau mean sea level
For i = 2 To 19 Step 2
j = i / 2
A = MatrixX(i): .Offset(j, 0) = A
B = MatrixX(i + 1): .Offset(j, 1) = B
ph = Atn(B / A) 'kwadran I
If A < 0 Then
ph = ph + pi 'kwadran II dan III
Else
If B < 0 Then ph = ph + 2 * pi 'kwadran IV
End If
'phase dikonversi ke derajat
Phase(j) = ph * 180 / pi: .Offset(j, 2) = Phase(j)
H(j) = Sqr(A * A + B * B): .Offset(j, 3) = H(j)
Next i
End With
'6. Membandingkan muka air hasil pengukuran dengan muka air hasil hitungan
Dim Ht As Double, SumHCos As Double
Dim Cetak() As Double
ReDim Preserve Cetak(1 To UBound(MatrixL), 1 To 5)
Zo = MatrixX(1)
For i = LBound(MatrixL) To UBound(MatrixL)
Cetak(i, 5) = CDbl(i)'mencetak nomer urut
Cetak(i, 1) = Range(addFirstDate) + (i - 1) / 24'mencetak hari dan jam
Cetak(i, 2) = MatrixL(i)'mencetak muka air pengukuran
SumHCos = 0
For j = 1 To 9
SumHCos = SumHCos + H(j) * Cos(w(j) * i + (Phase(j) * pi / 180))'rumus harmonic
Next j
Ht = Zo + SumHCos 'rumus(1)
Cetak(i, 3) = Ht'cetak muka air hasil hitungan
Cetak(i, 4) = Ht - MatrixL(i)'muka air hitungan - pengukuran
Next i
Range("A90:E785") = Cetak
End Sub
Dengan menggunakan data perbandingan antara muka air hasil pengukuran dengan perhitungan, maka dapat dibuat graphic sebagai berikut:
dari hasil hitungan didapat standard deviasi sebesar 0.11meter.
Prediksi pasang surut (pasut) menggunakan rumus Harmonic Analysis dengan metode least square untuk mencari 9 constituents utama pembangkit pasut dan 9 phase.
Dalam referensi tertulis,
"The ideal tide curve for any given port is represented as an average height Z0 plus a sum of terms ("constituents'') each of which is of the form f(t) = H cos(at + \phi). The time t is measured in hours, and f comes out in feet. The numbers H,a, \phi are the amplitude, the speed and the phase of the constituent."
Karena data hasil bacaan pasut yang tersedia adalah dalam meter dan dicatat tiap jam, maka nilai f tidak lagi dinyatakan dalam feet, melainkan dalam meter.
Rumus kemudian dikembangkan menjadi:
dalam format matrix dapat ditulis:
Jika MatrixX telah didapat, maka dengan menggunakan rumus(4) dan rumus(5),amplitudo dan phase tiap constituent dapat dihitung.
INPUT (DATA HASIL PENGUKURAN MUKA AIR)
Hasil pengukuran muka air didapat dari pengamatan selama 29 hari dengan pencatatan tiap jam, dimulai jam 0:00 sampai jam 23:00 tanggal 1 Maret 2008 sampai jam 23:00 tanggal 29 Maret 2008 yang ditabelkan sebagai berikut:
PROSES PEMBACAAN DATA DAN PERHITUNGAN
Proses pembacaan dan perhitungan dengan menggunakan program Visual Basic Application (macro) yang ada di microsoft excel 2003.
Khusus untuk prosedure atau program perkalian matirx dan inverse didapat dari http://www.alglib.net/
Pada file excel spread sheet saya, data pengukuran muka air pertama (1.850m) ada di cell $C$10 sedangkan data terakhir (2.110m) ada di $Z$38, sehingga range data pengukuran adalah di $C$10:$Z$38.
Pastikan bahwa di bagian Declaration di awal program diset :
Option Explicit 'berguna untuk mendeteksi definisi variable
Option Base 1 'hitungan matrix dimulai dari 1, kalo tidak diset, maka default index matrix adalah 0
Sub PrediksiPasutDenganLeastSquare()
'1. membaca data pengukuran muka air
dim cr as Range, rgData as Range
set rgData=Activesheet.Range("$C$10:$Z$38")
For Each cr In rgData.Columns(1).Cells
i = i + 1
j = 1 + (i - 1) * 24
ReDim Preserve MatrixL(1 To j + 23)
For jm = 0 To 23
MatrixL(j + jm) = cr.Offset(0, jm).Value
Next jm
Next
'2. membaca atau set variable periode tiap Constituents
Dim w(1 To 9) As Double 'periode dari 9 Constituents pasut
Dim pi As Double
pi = 4 * Atn(1)
w(1) = 2# * pi / 12.4206 'M2
w(2) = 2# * pi / 12# 'S2
w(3) = 2# * pi / 12.6582 'N2
w(4) = 2# * pi / 11.9673 'K2
w(5) = 2# * pi / 23.9346 'K1
w(6) = 2# * pi / 25.8194 'O1
w(7) = 2# * pi / 24.0658 'P1
w(8) = 2# * pi / 6.2103 'M4
w(9) = 2# * pi / 6.1033 'MS4
'3. membuat atau membetuk MatrixA() atau Matrix coeffisien
Dim it As Integer,MatrixA() as Double
ReDim Preserve matrixA(1 To UBound(MatrixL), 1 To 19)
For i = 1 To UBound(MatrixL)
matrixA(i, 1) = 1
it = i
For j = 1 To 9
matrixA(i, 2 * j) = Cos(w(j) * it)
matrixA(i, 2 * j + 1) = -Sin(w(j) * it)
Next j
Next i
'4. proses least square
Dim weight() As Double: ReDim Preserve weight(UBound(MatrixL, 1))
For i = LBound(weight) To UBound(weight): weight(i) = 1: Next i 'matrix bobot=matrix identitas
MatrixX() = clsLSQ.LSPAR(matrixA, MatrixL, weight)
'clsLSQ adalah class yang saya buat untuk proses perhitungan least square, tidak dibahas di session ini
'5. menghitung amplitudo dan phase 9 konstituent dan sekaligus menampilkan hasilnya di excel
Dim A As Double, B As Double, H(1 To 9) As Double, Phase(1 To 9) As Double
Dim ph As Double, Zo As Double
Const addPrint As String = "H66"
With Range(addPrint)
.Offset(0, 3) = MatrixX(1)'mencetak Zo atau mean sea level
For i = 2 To 19 Step 2
j = i / 2
A = MatrixX(i): .Offset(j, 0) = A
B = MatrixX(i + 1): .Offset(j, 1) = B
ph = Atn(B / A) 'kwadran I
If A < 0 Then
ph = ph + pi 'kwadran II dan III
Else
If B < 0 Then ph = ph + 2 * pi 'kwadran IV
End If
'phase dikonversi ke derajat
Phase(j) = ph * 180 / pi: .Offset(j, 2) = Phase(j)
H(j) = Sqr(A * A + B * B): .Offset(j, 3) = H(j)
Next i
End With
'6. Membandingkan muka air hasil pengukuran dengan muka air hasil hitungan
Dim Ht As Double, SumHCos As Double
Dim Cetak() As Double
ReDim Preserve Cetak(1 To UBound(MatrixL), 1 To 5)
Zo = MatrixX(1)
For i = LBound(MatrixL) To UBound(MatrixL)
Cetak(i, 5) = CDbl(i)'mencetak nomer urut
Cetak(i, 1) = Range(addFirstDate) + (i - 1) / 24'mencetak hari dan jam
Cetak(i, 2) = MatrixL(i)'mencetak muka air pengukuran
SumHCos = 0
For j = 1 To 9
SumHCos = SumHCos + H(j) * Cos(w(j) * i + (Phase(j) * pi / 180))'rumus harmonic
Next j
Ht = Zo + SumHCos 'rumus(1)
Cetak(i, 3) = Ht'cetak muka air hasil hitungan
Cetak(i, 4) = Ht - MatrixL(i)'muka air hitungan - pengukuran
Next i
Range("A90:E785") = Cetak
End Sub
Dengan menggunakan data perbandingan antara muka air hasil pengukuran dengan perhitungan, maka dapat dibuat graphic sebagai berikut:
dari hasil hitungan didapat standard deviasi sebesar 0.11meter.
