Сергеев В.Л. Интегрированные системы идентификации

Подождите немного. Документ загружается.

181

'1.4. Начальные значения вектора параметров AL

ALo(1) = 11

ALo(2) = 0.22

ALo(3) = 1.6

' 1.5. Задание переменных для адаптации ИСИ

kr = 0.0025 'Задание значения управляющего параметра

h = 1 'параметр шага в методе Гаусса–Ньютона

' Задание начальных значений параметров модели AL = ALo

AL(1) = ALo(1)

AL(2) = ALo(2)

AL(3) = ALo(3)

' 2. Алгоритм синтеза оптимальных оценок параметров КВД с использованием метода

' Гаусса–Ньютона

i2 = 0 '- переменная внешнего цикла

Do ' Условный оператор для уточнения вектора параметров AL(i), i = 1, 2, 3

' 2.1. Формирование матрицы частных производных – D(n,m) модели добычи нефти

For i = 1 To n

d(i, 1) = Exp(-AL(2) * i) * Exp(AL(3) * Log(i))

d(i, 2) = -AL(1) * i * Exp(-AL(2) * i) * Exp(AL(3) * Log(i))

d(i, 3) = AL(1) * AL(3) * Exp(-AL(2) * i) * Exp((AL(3) – 1) * Log(i))

Next i

' 2.2. Транспонирование матрицы частных производных – D(n,m)

For i = 1 To m

For j = 1 To n

DT(i, j) = d(j, i)

Next j

Next i

'2.3. Вычисление произведения матриц DTD

For i = 1 To m

For j = 1 To m

DTD(i, j) = 0

Next j

Next i

For i = 1 To m

For j = 1 To m

For k = 1 To n

DTD(i, j) = DTD(i, j) + DT(i, k) * d(k, j)

Next k

Next j

Next i

182

' 2.4. Формирование вектора-строки Ds

For i = 1 To m

DTs(1, i) = 0

Next i

For i = 1 To T

DTs(1, 1) = DTs(1, 1) + Exp(-AL(2) * i) * Exp(AL(3) * Log(i))

DTs(1, 2) = DTs(1, 2) + (-AL(1) * i * Exp(-AL(2) * i) * Exp(AL(3) * Log(i)))

DTs(1, 3) = DTs(1, 3) + AL(1) * AL(3) * Exp(-AL(2) * i) * Exp((AL(3) – 1) * Log(i))

Next i

'2.5. Транспонирование вектора Ds

For i = 1 To m

Ds(i, 1) = DTs(1, i)

Next i

' 2.6. Расчет матрицы частных производных DTsDs

For i = 1 To m

For j = 1 To m

DTsDs(i, j) = 0

Next j

Next i

For i = 1 To m

For j = 1 To m

DTsDs(i, j) = DTsDs(i, j) + Ds(i, 1) * DTs(1, j)

Next j

Next i

'2.7. Формирование матрицы планирования ИСИ A = DTD+ kr*DTsDs

For i = 1 To m

For j = 1 To m

A(i, j) = DTD(i, j) + kr * DTsDs(i, j)

Next j

Next i

'2.8. Формирование вектора-столбца невязок между фактическими и

' модельными значениями добычи нефти E

For i = 1 To n

F(i) = AL(1) * Exp(-AL(2) * i) * Exp(AL(3) * Log(i))

E(i, 1) = y(i) – F(i)

Next i

'2.9. Формирование вектора-столбца DTE

For i = 1 To m

For k = 1 To n

DTE(i, 1) = DT(i, k) * E(k, 1)

Next k

Next i

183

'2.10. Вычисление невязки между средней экспертной оценкой извлекаемых запасов

' и фактическими извлекаемыми запасами s – Es = MS_– s

s = 0

For i = 1 To T

s = s + AL(1) * Exp(-AL(2) * i) * Exp(AL(3) * Log(i))

Next i

Es = MS_ – s

'2.11. Вычисление вектора-столбца DTsEs

For i = 1 To m

DTsEs(i, 1) = DTs(1, i) * Es

Next i

'2.12. Вычисление вектора правых частей ИСМ( B = DTE+kDTsEs)

For i = 1 To m

B(i) = kr * DTsEs(i, 1) + DTE(i, 1)

Next i

' 2.13. Решение СЛУ по схеме Холецкого(процедура ELS_hol_1).

Call ELS_hol_1(B1, C, Y1, A, B, DeltaAL, m)

'Расчет функционала качества QLo на шаге L-1

i1 = 1

QLh(i1) = 0

For i = 1 To n

QLh(i1) = QLh(i1) + (y(i) – F(i)) * (y(i) – F(i))

Next i

' 2.14. Определение вектора параметров AL(L) = AL(l-1)+h*Delta AL(L-1)

For i = 1 To m

ALh(i) = AL(i)

Next i

' 2.15. Проверка условия сходимости метода Гаусса–Ньютона

For i = 1 To m

AL(i) = ALh(i) + h * DeltaAL(i)

Next i

'Расчет функционала качества QL1 на шаге L

QLh(i1 + 1) = 0

For i = 1 To n

F(i) = AL(1) * Exp(-AL(2) * i) * Exp(AL(3) * Log(i))

QLh(i1 + 1) = QLh(i1 + 1) + (y(i) – F(i)) * (y(i) – F(i))

Next i

' Дробление шага

i1 = i1 + 1

h = h / 2

' Проверка условия сходимости процедуры Гаусса–Ньютона

Loop While (QLh(i1 – 1) > = QLh(1)) And (i1 < = 10) i2 = i2 + 1

184

QL(i2 – 1) = QLh(1)

QL(i2) = QLh(i1 – 1)

' Проверка условия точности оценки параметров AL

Loop While ((QL(i2 – 1) – QL(i2)) / QL(i2) > = 0.1) And (i2 < = 3)

'3. Расчет точности оценок и вывод результатов

'3.1. Расчет оценки извлекаемых запасов при AL = AL*

SAL = 0

For i = 1 To T

SAL = SAL + AL(1) * Exp(-AL(2) * i) * Exp(AL(3) * Log(i))

Next i

' Расчет относительной ошибки оценки извлекаемых запасов в зависимости от числа лет разработки

DeltaSAL = Abs((SAL – S_) / S_)

' 3.2 Расчет прогноза добычи нефти на L лет

YL(n + L) = AL(1) * Exp(-AL(2) * (n + L)) * Exp(AL(3) * Log(n + L))

' Расчет относительной ошибки прогнозных значений добычи нефти

DeltaYL1 = ALy(1) * Exp(-ALy(2) * (n +L)) * Exp(ALy(3) * Log(n + L))

DeltaYL = Abs((YL(n + L) – DeltaYL1) /DeltaYL1)

Next n ‘ окончание цикла по числу лет истории разработки месторождения

End Sub

_______________________________________________________________________________

' Процедура решения СЛУ методом Холецкого

Private Sub ELS_hol_1(B1() As Single, C() As Single, Y1() As Single, A() As Single, B() _

As Single, DeltaAL() As Single, m As Integer)

For i = 1 To m

B1(i, i) = 1

C(i, 1) = A(i, 1)

Next i

For i = 2 To m

B1(1, i) = A(1, i) / A(1, 1)

Next i

' Формирование треугольных матриц С и В1

For j = 2 To m

' Формирование нижней треугольной матрицы С

For i = j To m sum = 0 For k = 1 To j – 1

sum = sum + C(i, k) * B1(k, j) Next k

C(i, j) = A(i, j) – sum

Next i

' Формирование верхней треугольной матрицы

В1

For i = j + 1 To m

sum = 0

For k = 1 To j – 1

sum = sum + C(j, k) * B1(k, i)

Next k

B1(j, i) = (A(j, i) – sum) / C(j, j)

Next i

185

Next j

'Прямой ход

Y1(1) = B(1) / C(1, 1)

For i = 1 To m sum = 0

For k = 1 To i – 1

sum = sum + C(i, k) * Y1(k)

Next k

Y1(i) = (B(i) – sum) / C(i, i)

Next i

' Обратный ход

DeltaAL(m) = Y1(m)

For i = 1 To m – 1 sum = 0

For k = m – i + 1 To m

sum = sum + B1(m – i, k) * DeltaAL(k)

Next k

DeltaAL(m – i) = Y1(m – i) – sum

Next i

2. Текст программы статистического моделирования и определения

параметров интегрированной системы модели кривой восстановления

забойного давления

Sub Example_2_В5a()

' Описание переменных динамических массивов данных

Dim n As Integer ' Число замеров КВД (nmax = Tmax*delta_t)

Dim m As Integer ' Число параметров модели КВД

Dim Tmax As Integer ' Максимальное число замеров КВД

Dim delta_t As Single ' Интервал между замерами КВД

Dim T1 As Single ' Максимальное время снятия КВД

Dim m1 As Integer ' Число экспертных оценок пластового давления

Dim m21 As Integer ' Число экспертных оценок первого параметра AL1 КВД

Dim m22 As Integer ' Число экспертных оценок второго параметров КВД

Dim m23 As Integer ' Число экспертных оценок третьего параметров КВД

Dim m3 As Integer ' Число экспертных оценок накопленной добычи жидкости

m = 3

L = 3

T1 = 150

delta_t = 1

Tmax = Int(T1 / delta_t) m1 = 5

m21 = 5

m22 = 5

m23 = 5

m3 = 2

For n = 1 To 15 ' Организация цикла по числу замеров КВД

'_______________________________________________________________________________

186

' Описание переменных и констант

ReDim y(n) As Single ' Имитируемый вектор измерений КВД

ReDim ALy(m) As Single ' Точные значения параметров модели КВД

ReDim AL(m) As Single ' Оценки параметров модели КВД

ReDim ALo(m) As Single ' Начальные значения параметров модели КВД

ReDim AL_(m) As Single ' Вектор средних значений экспертных оценок параметров КВД

ReDim AL_1(m21) As Single ' Вектор экспертных оценок первого параметра КВД

ReDim AL_2(m22) As Single ' Вектор экспертных оценок второго параметра КВД

ReDim AL_3(m22) As Single ' Вектор экспертных оценок третьего параметра КВД

ReDim ALpo(m) As Single

ReDim ALh(m) As Single ' Вспомогательный вектор параметров AL при дроблении шага h

' в методе Гаусса–Ньютона

ReDim ALp(m) As Single ' Вспомогательный вектор для преобразования параметров КВД

ReDim ND1(Tmax) As Single ' Вектор случайных величин для задания ошибок измерений КВД

ReDim ND2(m1) As Single ' Вектор случайных величин в модели пластового давления

ReDim ND21(m21) As Single ' Вектор случайных величин для задания ошибок параметра КВД AL(1)

ReDim ND22(m22) As Single ' Вектор случайных величин для задания ошибок параметра КВД AL(2)

ReDim ND23(m23) As Single ' Вектор случайных величин для задания ошибок параметра КВД AL(3)

ReDim ND3(m3) As Single ' Вектор случайных величин в модели накопленной жидкости в скважине

ReDim Pw_(m1) As Single ' Вектор экспертных оценок пластового давления

ReDim Sж_(m3) As Single ' Вектор экспертных оценок накопленной жидкости в скважине

ReDim F(n) As Single ' Вектор значений выхода модели КВД

ReDim Fo(n) As Single ' Вектор значений выхода модели КВД при AL = ALo

ReDim dFy(Tmax) As Single ' Производная модели КВД по времени

ReDim FL(Tmax) As Single ' Вектор значений логарифма в моменты времени i*delta_t,i = 1,2...

ReDim d(n, m) As Single ' Матрица частных производных (ЧП) по параметрам модели КВД

ReDim DT(m, n) As Single ' Транспонированная матрица ЧП

ReDim DTD(m, m) As Single ' Матрица планирования для модели КВД

ReDim dp(m, 1) As Single ' Вектор ЧП по параметрам AL в модели пластового давления

ReDim DTp(1, m) As Single ' Транспонированный вектор ЧП dp(m,1)

ReDim DTp1(Tmax, m) As Single

ReDim DTpDp(m, m) As Single ' Матрица частных производных модели априорной информации

' о пластовом давлении

ReDim Gamma(m) As Single ' Индикаторная матрица в модели экспертных оценок параметров КВД

ReDim G(m) As Single ' Индикаторная матрица в модели экспертных оценок накопленной

‘ жидкости

ReDim DS(m, 1) As Single ' Вектор ЧП по параметрам AL в модели накопленной жидкости в скважине

ReDim DTS(1, m) As Single ' Транспонированный вектор ЧП DS(m)

ReDim DSDTS(m, m) As Single ' Матрица ЧП по параметрам модели накопленной жидкости

ReDim Ep(n, 1) As Single ' Вектор столбец невязок по забойным давления ЕP = P*з-Pз

ReDim EAL_(m) As Single ' Вектор столбец AL_ – AL

ReDim DTE(m, 1) As Single ' Вектор правых частей системы линейных уравнений(СЛУ)

ReDim DTpEp(m, 1) As Single

ReDim A(m, m) As Single ' Матрица СЛУ для интегрированной системы моделей КВД

ReDim B(m) As Single ' Вектор правых частей СЛУ для ИСМ КВД

ReDim DeltaAL(m) As Single ' Вектор приращений параметров AL(m)

187

ReDim DeltaALy(m) As Single ' Относительная ошибка оценок параметров КВД

ReDim RAL(m) As Single ' Вспомогательный вектор для формирования априорной информации

' о параметрах КВД AL(m)

ReDim Rs(m) As Single ' Вспомогательный вектор для формирования априорной информации

' о накопленной добыче нефти S_

Dim P_ As Single ' Среднее значение экспертных оценок пластового давления

Dim S_ As Single ' Среднее значение экспертных оценок накопленной добыче нефти

Dim qo As Single ' Дебит скважины до ее остановки

Dim PAL As Single ' Оценка пластового давления

Dim P_ALy As Single ' Точное значение пластового давления

Dim P_dFy As Single ' Значение пластового давления через интеграл от производной КВД

Dim DeltaPAL As Single ' Относительная ошибка оценки пластового давления

Dim kr1 As Single ' Параметр регуляризации

Dim kr2 As Single ' Параметр регуляризации

Dim h As Single ' Параметр шага h в методе

Гаусса–Ньютона

Dim QLh(15) As Single ' Функционал качества (для проверки условия сходимости метода ГН)

Dim QL(15) As Single ' Функционал качества (для критерия точности оценки параметров)

Dim Ep_ As Single ' Невязка по пластовому давлению Ep_ = P_ – P(T,AL)

Dim Es As Single ' Невязка по накопленной жидкости Es = S_ – S(T,AL(2))

Dim DM(10) As Single ' Оценки детерминированных моментов КВД

Dim DMT(10) As Single ' Точные значения детерминированных моментов КВД

Dim ND(500) As Single ' Базовый вектор псевдослучайных чисел – N(0,1)

Dim i1 As Integer ' Число дробления параметра h

Dim i2 As Integer ' Переменная номера шага в мет

оде Гаусса–Ньютона

'_______________________________________________________________________________

ReDim B1(m, m) As Single ' В1,С,Y1 – массивы для решения СЛУ

ReDim C(m, m) As Single ' в процедуре ELS_hol_1

ReDim Y1(m) As Single ' C – нижняя треугольная матрица, B1 – верхняя треугольная матрица

'_______________________________________________________________________________

'1. Формирование интегрированной системы моделей КВД

'1.1. Формирование векторов псевдослучайных чисел

'1.1.1. Формирование базового вектора ND(500) псевдослучайных чисел – N(0,1)

Randomize (78) 'Начальные значения датчика равномерных псевдослучайных величин Rnd()

Rnd (-506)

For i = 1 To 500 s = 0

For j = 1 To 6 s = s + Rnd()

Next j

ND(i) = (s – 3) * 1.4 Worksheets(Лист 28).Cells(i, 32) = ND(i)

Next i

'1.1.2. Формирование вектора случайных величин для имитации ошибок измерений КВД

For i = 1 To Tmax

ND1(i) = ND(i)

Next i

188

'1.1.3. Формирование вектора случайных величин для имитации экспертных оценок

'пластового давления

For i = 1 To m1

ND2(i) = ND(i + Tmax)

Next i

'1.1.4. Формирование вектора случайных величин для имитации экспертных оценок

' параметра AL(1) модели КВД

For i = 1 To m21

ND21(i) = ND(i + Tmax + m1)

Next i

'1.1.5. Формирование вектора случайных величин для имитации экспертных оценок

' параметра AL(2) модели КВД

For i = 1 To m22

ND22(i) = ND(i + Tmax + m1 + m22)

Next i

'1.1.6. Формирование вектора случайных величин для имитации экспертных оценок

' параметра AL(3) модели КВД

For i = 1 To m23

ND23(i) = ND(i + Tmax + m1 + m22 + m23)

Next i

'1.1.7. Формирование вектора случайных величин для имитации экспертных оценок

' накопленной добычи жидкости в скважине после ее остановки

For i = 1 To m3

ND3(i) = ND(i + Tmax + m1 + m22 + m23 + m3)

Next i

'1.2. Формирование вектора измерений КВД Y*

' Задание точных значений параметров КВД

ALy(1) = 12

ALy(2) = 0.7

ALy(3) = 2

qo = 1 ' Задание значения дебита до остановки скважины

ALp(3) = Log(ALy(3))

Worksheets(Лист 28).Cells(13, 12) = ALy(1)

Worksheets(Лист 28).Cells(13, 13) = ALy(2)

Worksheets(Лист 28).Cells(13, 14) = ALp(3)

' Расчет значений логарифма FL = Lnt в точках ti,i = 1,2,...,Tmax

For i = 1 To Tmax

FL(i) = Log(i * delta_t)

Next i

'Расчет вектора измерений КВД Y*

yo = 180 ' Начальное значение забойного давления Рз(to)

с1 = 0.005 ' Относительный уровень ошибок измерений забойного давления

For i = 1 To n

y(i) = yo + ALy(1) * qo * (1 – Exp(-ALy(2) * i * delta_t)) * (ALp(3) + FL(i))

y(i) = y(i) + с1 * y(i) * ND(i)

Worksheets(Лист 28).Cells(i + 1, 3) = y(i)

Next i

189

' 1.3. Формирование экспертных оценок пластового давления

'1.3.1. Приближенный расчет пластового давления по модели S_dFy

P_dFy = yo ' начальное значение забойного давления

For i = 1 To Tmax

P_dFy = P_dFy + ALy(1) * ALy(2) * qo * Exp(-ALy(2) * i * delta_t) * (Log(ALy(3)) + FL(i)) _

+ ALy(1) * qo * (1 – Exp(-ALy(2) * i * delta_t)) / (i * delta_t)

Next i

Worksheets(Лист 28).Cells(9, 17) = P_dFy

' 1.3.2. Точное значение пластового давления при t = Tmax (248,45 ат)

P_ALy = yo + ALy(1) * qo * (1 – Exp(-ALy(2) * Tmax)) * (ALp(3) + FL(Tmax))

Worksheets(Лист 28).Cells(11, 17) = P_ALy

' 1.3.3. Моделирование экспертных оценок пластового давления

c2 = 0.02 'Относительный уровень ошибок экспертных оценок пластового давления – 5 %

P_ = 0

For i = 1 To m1

Pw_(i) = P_ALy + c2 * P_ALy * ND2(i)

P_ = P_ + Pw_(i)

Worksheets(Лист 28).Cells(i + 19, 20) = Pw_(i)

Next i

P_ = P_ / m1

Worksheets(Лист 28).Cells(20, 21) = P_

Worksheets(Лист 28).Cells(20, 22) = c2

'1.4. Формирование экспертных оценок о параметрах КВД

'1.4.1. Формирование экспертных оценок о параметре AL(1)

c21 = 0.2 'Относительный уровень ошибок экспертных оценок пластового давления – 5 %

AL_(1) = 0

For i = 1 To m21

AL_1(i) = ALy(1) + c21 * ALy(1) * ND21(i)

AL_(1) = AL_(1) + AL_1(i)

Worksheets(Лист 28).Cells(i + 42, 14) = AL_1(i)

Next i

AL_(1) = AL_(1) / m21

Worksheets(Лист 28).Cells(43, 17) = AL_(1)

Worksheets(Лист 28).Cells(43, 20) = c21

'1.4.2. Формирование экспертных оценок о параметре AL(2)

c22 = 0.2 'Относительный уровень ошибок экспертных оценок пластового давления – 5 %

AL_(2) = 0

For i = 1 To m22

AL_2(i) = ALy(2) + c22 * ALy(2) * ND22(i)

AL_(2) = AL_(2) + AL_2(i)

Worksheets(Лист 28).Cells(i + 42, 15) = AL_2(i)

Next i

AL_(2) = AL_(2) / m22

Worksheets(Лист 28).Cells(43, 18) = AL_(2)

Worksheets(Лист 28).Cells(43, 21) = c21

190

'1.4.3. Формирование экспертных оценок о параметре AL(3)

c23 = 0.2 'Относительный уровень ошибок экспертных оценок пластового давления – 5 %

AL_(3) = 0

For i = 1 To m23

AL_3(i) = ALy(3) + c23 * ALy(3) * ND23(i)

AL_(3) = AL_(3) + AL_3(i)

Worksheets(Лист 28).Cells(i + 42, 16) = AL_3(i)

Next i

AL_(3) = AL_(3) / m23

Worksheets(Лист 28).Cells(43, 19) = AL_(3)

Worksheets(Лист 28).Cells(43, 22) = c23

'1.4.4. Формирование экспертных оценок о накопленной жидкости в скважине S_

' S0_ = qo/ALy(2) – Точное значение накопленной жидкости в скважине

S0_ = qo / ALy(2)

Worksheets(Лист 28).Cells(32, 23) = S0_

c3 = 0.2 'Относительный уровень ошибок экспертных оценок накопленной жидкости – 20 %

S_ = 0

For i = 1 To m3

Sж_(i) = S0_ + c3 * S0_ * ND3(i)

S_ = S_ + Sж_(i)

Worksheets(Лист 28).Cells(i + 31, 20) = AL_3(i)

Next i

S_ = S_ / m3

Worksheets(Лист 28).Cells(32, 21) = S_

Worksheets(Лист 28).Cells(32, 22) = c3

'2. Алгоритм синтеза оценок параметров КВД с использованием метода Гаусса–Ньютона

' Начальные значения вектора параметров AL

ALo(1) = 13

ALo(2) = 0.8

ALo(3) = 2.2

ALpo(3) = Log(ALo(3))

Worksheets(Лист 28).Cells(13, 8) = ALo(1)

Worksheets(Лист 28).Cells(13, 9) = ALo(2)

Worksheets(Лист 28).Cells(13, 10) = ALpo(3)

' Расчет вектора модели КВД при AL = ALo

For i = 1 To n

Fo(i) = yo + ALo(1) * qo * (1 – Exp(-ALo(2) * i * delta_t)) * (ALpo(3) + FL(i))

Worksheets(Лист 28).Cells(i + 1, 2) = Fo(i)

Next i

' Расчет пластового давления при начальных значениях параметров AL =

ALo

P_ALo = yo + ALo(1) * qo * (1 – Exp(-ALo(2) * Tmax)) * (ALp(3) + FL(Tmax))

Worksheets(Лист 28).Cells(14, 17) = P_ALo

'Задание начальных значений параметров модели AL = ALo

AL(1) = ALo(1)

AL(2) = ALo(2)