Приложение b. Доказательство
Разложим
нелинейную функцию
в ряд Тейлора в окресности некоторой
точки
(начальное приближение параметров),
ограничиваясь членами первого порядка
малости:
С
учетом этого разложения критерий
качества модели исследуемого объекта
представляется
в виде
где
– матрица частных производных по
параметрам
.
Переходя от нелинейной интегрированной системы модолей к линейной интегрированной системк моделей
относительно
приращений вектора параметров
;
.
Линеаризованной интегрированной системе моделей соответствует комбинированный критерий качества
Определим
приращение вектора параметра
относительно начального приближения
и
получим последующие приближения
параметров по схеме Ньютона.
Необходимо
вычислить производные от функционала
по параметрам
и приравнять к нулю. В результате получим
систему линейных уравнений
Приложение с. Листинг программы
Таблица 1 – Перечень процедур и функций модуля 1.
|
Sub Example() |
Процедура, в которой производится считывание данных и вызов всех ниже перечисленных функций для рассчета значений КВД |
|
Function Func(ByVal t As Single, A() As Single) As Double |
Функция, задающая модель построения КВД |
|
Private Sub SimplexMethod(AL() As Single, ByVal n As Single, y() As Single, ByVal m As Integer, kr As Single, MS_ As Single, ExtractionTime As Integer) |
Процедура расчета параметров методом Деформированного многогранника |
|
Public Function ArgMinOpt(AL() As Single, ByVal n As Integer, y() As Single) As Single |
Функция расчета показателя качества управляющего параметра |
|
Function ArgMin(AL() As Single, ByVal n As Integer, y() As Single, kr As Single, MS_ As Single, t As Integer) As Single |
Комбинированный показатель качества ИСМ добычи нефти |
Таблица 2 – Перечень процедур и функций модуля №2.
|
Sub Example_2() |
Процедура, в которой производится считывание данных и вызов всех ниже перечисленных функций для рассчета значений КВД |
|
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) |
Процедура решения СЛУ методом Холецкого |
|
Function Func1(ByVal t As Single, A() As Single) As Double |
Функция, задающая модель построения КВД |
|
Public Sub GN(B1() As Single, C() As Single, Y1() As Single, A() As Single, B() _ As Single, DeltaAL() As Single, m As Integer, ByVal n As Integer, y() As Single, F() As Single, AL() As Single, _ ALh() As Single, kr As Single, MS_ As Single, ByVal t As Integer) |
Процедура определения параметров ИСМ добычи нефти методом Гаусса-Ньютона |
|
Function ArgMin(AL() As Single, ByVal n As Integer, y() As Single, kr As Single, t As Integer) As Single |
Вычисление функционала качества |
|
Option Explicit Sub Example() 'программа прогноза добычи нефти с учетом априорной информации и экспертных оценок об извлекаемых запасах 'используются методы оптимизации: Нелдера и Мида (деформируемого многогранника) и Гаусса-Ньютона 'модель добычи нефти имеет вид: y(e)= A(1) * Exp(-A(2) * t) * t ^ A(3). ' Dim n As Integer 'время разработки (колличество интервалов) Dim m As Integer 'число параметров модели объекта Dim L As Integer 'переменная времени прогноза добычи нефти For n = 1 To 21 'организация цикла по количеству интервалов времении разработки месторождения m = 3 L = 3
'Описание переменных
ReDim y(n + L) As Single 'Имитируемый вектор фактической добычи нефти ReDim YL(n + L) As Single 'Вектор прогноза добычи нефти на L моментов времени ReDim ALy(m) As Single 'Точные значения параметров модели добычи нефти ReDim AL(m) As Single 'Оценки параметров модели добычи нефти ReDim ALo(m) As Single 'Начальные значения параметров модели добычи нефти ReDim Fo(n) As Single 'Вектор значений добычи нефти при начальных значениях параметров модели AL=ALo ReDim F(n) As Single 'Вектор модельных значений добычи нефти
Dim S_ As Single 'Фактические извлекаемые запасы нефти Dim MS_ As Single 'Экспертная оценка извлекаемых запасов Dim t As Integer 'время разработки Dim SAL As Single 'Оценка извлекаемых запасов Dim DeltaSAL As Single 'Относительная ошибка оценки извлекаемых запасов Dim kr As Single 'управляющий параметр (параметр регуляризации) Dim c1 As Single 'Уровень ошибок измерений добычи нефти Dim c2 As Single 'Уроввень ошибок экспертных оценок извлекаемых запасов нефти
Dim i As Integer Dim j As Integer
Dim DeltaYL As Single Dim DeltaYL1 As Single Dim DeltaMS_ As Single Dim temp As Single
'1 ФОРМИРОВАНИЕ ИСМ ДОБЫЧИ НЕФТИ '1.1 Формирование вектора фактических значений добычи нефти y*(t) t = 30 'время разработки месторождения
'Ввод данных КВД скв 121 For i = 1 To n y(i) = Worksheets("Лист1").Cells(i + 1, 4) Next i
'1.2 Формирование экспертных оценок извлекаемых запасов '1.2.1 расчет и вывод оценки начальной емкости рынка S S_ = 180 Worksheets("Лист1").Cells(34, 5) = S_
' 1.2.2 Экспертная оценка начальной емкости рынка с относительной ошибкой с2 c2 = 0.03 MS_ = S_ + S_ * c2 Worksheets("Лист1").Cells(5 + n, 5) = MS_
'1.3 Расчет вектора значений добычи нефти на основе модели при AL=ALo '1.3.1 начальные значения вектора параметров AL ALo(1) = 155 ALo(2) = 0.45 ALo(3) = 0.1 For i = 1 To n Fo(i) = Func(i, ALo) Worksheets("Лист1").Cells(i + 1, 2) = Fo(i) 'Вывод вектора параметров на лист Next i
'1.3.2 вывод начальных значений параметров модели добычи нефти Worksheets("Лист1").Cells(13, 8) = ALo(1) Worksheets("Лист1").Cells(13, 9) = ALo(2) Worksheets("Лист1").Cells(13, 10) = ALo(3) 'задание начальных значений параметров модели AL=ALo AL(1) = ALo(1) AL(2) = ALo(2) AL(3) = ALo(3)
'2. ОПРЕДЕЛЕНИЕ ПАРАМЕТРОВ AL МОДЕЛИ ДОБЫЧИ НЕФТИ МЕТОДОМ ДЕФОРМИРОВАННОГО МНОГОГРАННИКА Dim LeftEdge As Single Dim RightEdge As Single Dim Center As Single Dim LeftVal As Single Dim RightVal As Single Dim CenterVal As Single Dim Iterations As Integer
' задание границ расчета параметров Iterations = 0 LeftEdge = 0.0001 RightEdge = 0.5
Call SimplexMethod(AL, n, y, m, LeftEdge, MS_, t) 'вызов функции по расчету параметров МДМ LeftVal = ArgMinOpt(AL, n, y) Call SimplexMethod(AL, n, y, m, RightEdge, MS_, t) RightVal = ArgMinOpt(AL, n, y) Do Center = (RightEdge + LeftEdge) / 2 Call SimplexMethod(AL, n, y, m, Center, MS_, t) CenterVal = ArgMinOpt(AL, n, y) If LeftVal > RightVal Then LeftEdge = Center LeftVal = CenterVal Else RightEdge = Center RightVal = CenterVal End If Iterations = Iterations + 1 Loop While (RightEdge - LeftEdge > 0.1) kr = (RightEdge + LeftEdge) / 2
'Расчет параметров модели добычи нефти при нулевых значения правляющего параметра kr Call SimplexMethod(AL, n, y, m, 0, MS_, t)
'Вывод значений управляющего параметра Worksheets("Лист1").Cells(37, 24) = "Beta1" Worksheets("Лист1").Cells(37 + n, 24) = kr Worksheets("Лист1").Cells(34 + n, 4) = kr
'Вывод значений добычи нефти полученных на основе модели при AL=AL* For i = 1 To n F(i) = Func(i, AL) Next i
'3. РАСЧЕТ ТОЧНОСТИ ОЦЕНОК И ВЫВОДА РЕЗУЛЬТАТОВ '3.1. Расчет оценки извлекаемых запасов при AL=AL* '3.1.1 вывод значений добычи нефти полученные на основе адаптированной модели For i = 1 To n Worksheets("Лист1").Cells(i + 1, 1) = F(i) Next i '3.1.2 вывод оптимальных оценок параметров модели добычи нефти Worksheets("Лист1").Cells(16, 8) = AL(1) Worksheets("Лист1").Cells(16, 9) = AL(2) Worksheets("Лист1").Cells(16, 10) = AL(3)
Worksheets("Лист1").Cells(12, 13) = AL(1) Worksheets("Лист1").Cells(12, 14) = AL(2) Worksheets("Лист1").Cells(12, 15) = AL(3)
'3.1.3 вывод оценки извлекаемых запасов SAL = Func(t, AL) Worksheets("Лист1").Cells(2, 5) = SAL
'3.1.4 расчет относительной ошибки оценки извлекаемых запасов в зависимости от числа лет разработки DeltaSAL = Abs((SAL - S_) / S_) Worksheets("Лист1").Cells(34 + n, 1) = SAL Worksheets("Лист1").Cells(34 + n, 2) = DeltaSAL Worksheets("Лист1").Cells(34 + n, 3) = n
' Расчет гидропроводности---------------------------------------------------------------------------- Dim q0 As Single Dim mu As Single Dim sig As Single
q0 = 124 * 11.57 * 1.16 * 0.86 ' дебит скважины до остановки(в пластовых условиях в см.куб/сек) mu = 2.6 ' вязкость нефти в пластовых условиях -сП) sig = q0 / (4 * 3.14 * AL(2) * 0.86) ' оценка гидропроводности пласта -Д*см./ сП Worksheets("Лист1").Cells(70 + n, 5) = sig
' Расчет пъезопроводности - Xn Dim bn As Single Dim mn As Single Dim hp As Single Dim bc As Single Dim Xn As Single Dim Pi As Single
bn = 0.00011 ' Коэффициент сжимаемости формации нефти bc = 0.00001 ' Коэффициент сжимаемости системы mn = 0.2 ' Пористость hp = 1000 ' Толщина пласта в сантиметрах Pi = 3.14 Xn = q0 / (4 * Pi * AL(1) * hp * (mn * bn + bc)) Worksheets("Лист1").Cells(70 + n, 6) = Xn
Dim Tmax As Single Dim nk As Single nk = 20 ' Номер конечной точки полной КВД Dim nk_ As Single nk_ = 20 ' Номер конечной точки полной КВД ReDim Pl(n) As Single ' Вектор оценок пластового давления ReDim P3m(nk) As Single ' Вектор значений забойного давления P3 полученный на основе модели КВД Dim Pl_exp As Single ' Экспертное значение пластового давления ReDim t1(nk) As Single ' Время снятия КВД T=log(t)
For i = 1 To nk t1(i) = Worksheets("Лист3").Cells(2 + i, 2) Next i
Pl_exp = Worksheets("Лист3").Cells(2, 20)
' Оценка пластового давления Tmax = 50000 Pl(n) = Pl_exp + AL(1) * (1 - Exp((-AL(2)) * Tmax)) ^ AL(3) Worksheets("Лист1").Cells(70 + n, 7) = Pl(n)
' Расчет модельных значений КВД For i = 1 To nk P3m(i) = AL(1) * (1 - Exp((-AL(2)) * t1(i))) ^ AL(3) Worksheets("Лист1").Cells(72 + i, 8) = P3m(i) Next i
Dim T_exp1 As Single ' Экспертное время восстановления КВД Dim tau As Single '
T_exp1 = Worksheets("Лист3").Cells(3, 21) tau = (T_exp1 - t1(nk)) / nk_
ReDim P3m_1(nk_) As Single
' Расчет модельных значений КВД (Прогнозирование) For i = 1 To nk_ P3m_1(i) = AL(1) * (1 - Exp((-AL(2)) * (t1(nk) + tau * i))) ^ AL(3) Worksheets("Лист1").Cells(92 + i, 8) = P3m_1(i) Next i
ReDim Pl_(n) As Single
' Оценка пластового давления (Прогнозирование) Pl_(n) = Pl_exp + AL(1) * (1 - Exp((-AL(2)) * T_exp1)) ^ AL(3) Worksheets("Лист1").Cells(90 + n, 7) = Pl_(n)
'====== Расчет детерминированных моментов- m0,m1,m2 (высисление интегралов методом прямоугольников) ==========
Dim m0 As Single Dim m0_ As Single Dim m1 As Single Dim m1_ As Single Dim m2 As Single Dim m2_ As Single Dim m0__ As Single Dim m1__ As Single Dim m2__ As Single Dim T_exp_1 As Single Dim tau_ As Single ReDim t1_(nk + 1) As Single ' Время снятия КВД в секундах T1 = t
For i = 1 To nk + 1 t1_(i) = Worksheets("Лист3").Cells(1 + i, 3) Next i
m0 = 0 For i = 1 To nk - 1 m0 = m0 + (Pl(n) - P3m(i)) * (t1_(i + 1) - t1_(i)) Next i
T_exp_1 = Worksheets("Лист3").Cells(2, 21) tau_ = (T_exp_1 - t1_(nk + 1)) / nk_
m0_ = 0 For i = 1 To nk_ m0_ = m0_ + (Pl_(n) - P3m_1(i)) * tau_ Next i
m1 = 0 For i = 1 To nk - 1 m1 = m1 + (Pl(n) - P3m(i)) * (t1_(i + 1) - t1_(i)) * t1_(i) Next i
m1_ = 0 For i = 1 To nk_ m1_ = m1_ + (Pl_(n) - P3m_1(i)) * tau_ * (t1_(nk) + i * tau_) Next i
m2 = 0 For i = 1 To nk - 1 m2 = m2 + (Pl(n) - P3m(i)) * (t1_(i + 1) - t1_(i)) * t1_(i) ^ 2 Next i
m2_ = 0 For i = 1 To nk_ m2_ = m2_ + (Pl_(n) - P3m_1(i)) * tau_ * (t1_(nk) + i * tau_) ^ 2 Next i
For i = 1 To nk_ m0__ = m0 + m0_ m1__ = m1 + m1_ m2__ = m2 + m2_ Next i
Worksheets("Лист1").Cells(70 + n, 20) = m0__ Worksheets("Лист1").Cells(70 + n, 21) = m1__ Worksheets("Лист1").Cells(70 + n, 22) = m2__
' Расчет гидропроводности sigm диагностическго коэффициента d методом детерминированных моментов и пъезопроводности Dim sigm__ As Single Dim pezm__ As Single Dim R As Single Dim d__ As Single R = 10 sigm__ = (4 / (5 * 3.14)) * q0 * m1__ / m0__ ^ 2 Worksheets("Лист1").Cells(70 + n, 24) = sigm__ pezm__ = ((6.4 * m0__ * R * R) / m1__) Worksheets("Лист1").Cells(70 + n, 25) = pezm__ d__ = m0__ * m2__ / m1__ ^ 2 Worksheets("Лист1").Cells(70 + n, 26) = d__
'---------------------------------------------------------------------------------------------------- Next n
End Sub
'Модель вычисления КВД Function Func(ByVal t As Single, A() As Single) As Double Func = A(1) * ((1 - Exp((-A(2)) * t)) ^ A(3)) End Function
'МЕТОД ДЕФОРМИРОВАННОГО МНОГОГРАННИКА Private Sub SimplexMethod(AL() As Single, ByVal n As Single, y() As Single, ByVal m As Integer, kr As Single, MS_ As Single, ExtractionTime As Integer) ReDim P(m + 1, m) As Single 'Точки симплекса ReDim tempP(m) As Single ReDim F(m + 1) As Double 'значения функционала в точках симплекса ReDim NewP(m) As Single 'новая вершина симплекса ReDim Pmin(m) As Single ReDim Ua(m) As Double 'точка - центр тяжести грани, относительно которой происходит отражение
ReDim L(m) As Double 'Вектор, характеризующий растояние от вершиы iMax до противоположной грани ReDim L0(m) As Double ' Dim NewF As Double 'знчение функционала в новой вершине Dim Fmin As Double Dim SimplexLength As Double 'длина симплекса Dim i As Integer Dim j As Integer Dim iMax As Integer 'индекс вершины с максимальным значением функционала Dim iMin As Integer 'индекс вершины с минимальным значением функционала Dim h As Single
Dim k As Single Const E = 0.01
'Создаем начальный симплекс h = 2 For i = 1 To m + 1 For j = 1 To m If (i = 1) Or ((i - 1) = j) Then P(i, j) = AL(j) Else P(i, j) = AL(j) + h End If Next j Next i
For i = 1 To m + 1 For j = 1 To m tempP(j) = P(i, j) Next j F(i) = ArgMin(tempP, n, y, kr, MS_, ExtractionTime) Next i
Do 'Ищем вершину с максимальным значением фунционала iMax = 1 For i = 1 To m + 1 If (F(iMax) < F(i)) Then iMax = i End If Next i 'Нахождение вектора расстояния от вершины с максимальным значением функционала до противоположной грани For j = 1 To m Ua(j) = 0 For i = 1 To m + 1 If (i <> iMax) Then Ua(j) = Ua(j) + P(i, j) End If Next i Ua(j) = Ua(j) / m L0(j) = Ua(j) - P(iMax, j) Next j k = 0.01 Fmin = 65536 'прямой поиск минимума вдоль L Do For j = 1 To m L(j) = L0(j) * k Next j For j = 1 To m NewP(j) = Ua(j) + L(j) If (NewP(j) < 0) Then NewP(j) = 0 End If Next j NewF = ArgMin(NewP, n, y, kr, MS_, ExtractionTime) If (NewF < Fmin) Then Fmin = NewF For j = 1 To m Pmin(j) = NewP(j) Next j End If k = k * 1.05 Loop While (k < 6)
SimplexLength = 0 For j = 1 To m SimplexLength = SimplexLength + L(j) * L(j) Next j SimplexLength = Sqr(SimplexLength) If F(iMax) > Fmin Then F(iMax) = Fmin For j = 1 To m P(iMax, j) = Pmin(j) Next j Else Exit Do End If
Loop While ((SimplexLength > E) And (Fmin > E)) For j = 1 To m AL(j) = Pmin(j) Next j End Sub
'Показатель качества уравляющего параметра kr Public Function ArgMinOpt(AL() As Single, ByVal n As Integer, y() As Single) As Single Dim i As Integer Dim S As Single Dim temp As Single S = 0 For i = 1 To n temp = (y(i) - Func(i, AL)) / y(i) S = S + temp * temp Next i ArgMinOpt = S End Function
'Комбинированный показатель качества ИСМ добычи нефти ArgMin = S + kr * S1 Function ArgMin(AL() As Single, ByVal n As Integer, y() As Single, kr As Single, MS_ As Single, t As Integer) As Single Dim i As Integer Dim S As Single Dim S1 As Single Dim temp As Single Dim SAL As Single
'Частный показатель качества модели добычи нефти S = 0 For i = 1 To n temp = (y(i) - Func(i, AL)) / y(i) S = S + temp * temp Next i SAL = Func(t, AL) 'Частный показатель качества модели извлекаемых запасов S1 = (SAL - MS_) * (SAL - MS_) / MS_ ^ 2 'Комбинированный показатель качества ArgMin = S + kr * S1 End Function |
|
Sub Example_2() 'Программа прогноза добычи нефти ' с учетом априорной информации и экспертных оценок об извлекаемых запасах ' с использованием метода оптимизации "Гаусса-Ньютона" ' и модели добычи нефти вида y(e)= A(1) * Exp(-A(2) * t) ^ A(3).. ' Вывод результатов на лист 2, 1
'Описание переменных динамических массивов данных Dim n As Integer ' Число интервалов разработки Dim m As Integer ' Число параметров модели объекта Dim L As Integer ' Переменная числа лет прогноза добычи нефти For n = 1 To 21 ' Организация цикла по числу лет истории разработки месторождения ' до максимального значения n=15 m = 3 ' Задание числа параметров модели объекта
'_______________________________________________________________________________________________ ' Описание переменных задачи
ReDim y(n) As Single ' Вектор фактических значений добичи нефти ReDim YL(n) As Single ' Вектор прогноза добычи нефти на L лет ReDim ALy(m) As Single ' Точные значения параметров модели объекта ReDim AL(m) As Single ' Оценки параметров модели объекта ReDim ALo(m) As Single ' Начальные значения параметров модели объекта ReDim ALh(m) As Single ' Вспомогательный вектор параметров AL при дроблении шага h ' в методе Гаусса-Нъютона ReDim F(n) As Single ' Вектор значений выхода модели объета ReDim Fo(n) As Single ' Вектор значений выхода модели объета при AL = ALo ReDim A(m, m) As Single ' Матрица СЛУ для интегрированной системы ' моделей(ИСМ) добычи нефти ReDim B(m) As Single ' Вектор правых частей СЛУ для ИСМ ДН ReDim DeltaAL(m) As Single ' Вектор прирощений параметров AL(m) Dim S_ As Single ' Извлекаемые запасы Dim MS_ As Single ' Экспертная оценка извлекаемых запасов Dim t As Integer ' число лет разработки Dim SAL As Single ' Оценка извлекаемых запасов Dim DeltaSAL As Single ' Относительная ошибка оценки извлекаемых запасов Dim kr As Single ' Управляющий параметр- параметр регуляризации Dim h As Single ' Параметр шага h в методе Гаусса-Нъютона Dim ND(51) As Single ' Вектор случайных чисел распределенных по нормальному закону N(0,1)
Dim DeltaYL1 As Single ' Вспомогательная переменная
Dim q0 As Single ' Вязкость нефти
' Массивы для решения системы линейных уравнений(СЛУ) __________________________________
ReDim B1(m, m) As Single ' В1,С,Y1 - массивы для решения СЛУ ReDim C(m, m) As Single ' в процедуре ELS_hol_1 ReDim Y1(m) As Single ' C - нижняя треугольная матрица,B1-верхняя треугольная матрица
'___________________________________________________________________________________________ ' Получения массича ND(i),i=1,151 псевдочслучайных чисел распределенных по нормальному закону N(0,1)
Randomize (78) 'Начальные значения датчика равномерных псевдослучайных величин Rnd() Rnd (-506)
For i = 1 To 51 S = 0 For j = 1 To 6 S = S + Rnd() Next j ND(i) = (S - 3) * 1.4 Next i
' 1. Формирование интегрированной системы моделей(ИСМ) добычи нефти____________________
' 1.1 Формирование вектора фактических значений добычи нефти
'Ввод данных КВД скв 121 For i = 1 To n y(i) = Worksheets("Лист1").Cells(i + 1, 4) Next i ALy(1) = 180 'параметры модели добычи нефти ALy(2) = 0.19 ALy(3) = 0.7
ALo(1) = 180 ALo(2) = 0.18 ALo(3) = 0.695 For i = 1 To n Fo(i) = Func1(i, ALo) Worksheets("Лист2").Cells(i + 1, 2) = Fo(i) Next i
' Параметры модели добычи нефти Worksheets("Лист2").Cells(13, 13) = ALy(1) Worksheets("Лист2").Cells(13, 14) = ALy(2) Worksheets("Лист2").Cells(13, 15) = ALy(3)
c1 = 0.05 ' Уровень "вариаци" фактических значений добычи нефти относительно ее модели
Worksheets("Лист2").Cells(13, 8) = ALo(1) Worksheets("Лист2").Cells(13, 9) = ALo(2) Worksheets("Лист2").Cells(13, 10) = ALo(3)
t = 50 'время разработки месторождения
' 1.2 Расчет и вывод оценки начальной емкости рынка S_
S_ = Func1(i, ALy) Worksheets("Лист2").Cells(35, 5) = S_
' 1.2. Экспертная оценка начальной емкости рынка с относительной ошибкой c2 c2 = 0.03 MS_ = S_ + S_ * c2 Worksheets("Лист2").Cells(11, 17) = MS_
h = 1 'параметр шага в методе Гаусса-Ньютона
'Задание начальных значений параметров модели AL=ALo AL(1) = ALo(1) AL(2) = ALo(2) AL(3) = ALo(3)
'2. Определение параметров AL модели добычи нефти методом Гаусса- Ньютона ' и управляююцего параметра методом дихотомии
Dim LeftEdge As Single Dim RightEdge As Single Dim Center As Single Dim LeftVal As Single Dim RightVal As Single Dim CenterVal As Single Dim Iterations As Integer
Iterations = 0 LeftEdge = 0.001 RightEdge = 0.5
' 2.1 Вызов процедуры Гаусса- Ньютона, формировыание матриц и решение СЛУ Call GN(B1, C, Y1, A, B, DeltaAL, m, n, y, F, AL, ALh, LeftEdge, MS_, t)
' Метод дихотомии LeftVal = ArgMin(AL, n, y, LeftEdge, t) Call GN(B1, C, Y1, A, B, DeltaAL, m, n, y, F, AL, ALh, RightEdge, MS_, t) RightVal = ArgMin(AL, n, y, RightEdge, t) Do Center = (RightEdge + LeftEdge) / 2 Call GN(B1, C, Y1, A, B, DeltaAL, m, n, y, F, AL, ALh, Center, S_, t) CenterVal = ArgMin(AL, n, y, Center, t) If LeftVal > RightVal Then LeftEdge = Center LeftVal = CenterVal Else RightEdge = Center RightVal = CenterVal End If Iterations = Iterations + 1 Loop While (RightEdge - LeftEdge > 0.1) kr = (RightEdge + LeftEdge) / 2
Call GN(B1, C, Y1, A, B, DeltaAL, m, n, y, F, AL, ALh, kr, MS_, t)
' Расчет параметров модели добычи нефти при нулевых значениях правляющего параметра kr ' Call GN(B1, C, Y1, A, B, DeltaAL, m, n, y, F, AL, ALh, 0, MS_, t)
' Вывод значений управляющего параметра Worksheets("Лист2").Cells(37, 24) = "Упр. парам." Worksheets("Лист2").Cells(37 + n, 24) = kr Worksheets("Лист2").Cells(34 + n, 4) = kr
'3. Расчет точности оценок и вывод рультатов '3.1. Расчет оценки извлекаемых запасов при AL=AL*
'Расчет оценки извлекаемых запасов при AL=AL* SAL = 0 ' For i = 1 To t SAL = SAL + Func1(t, AL) ' Next i Worksheets("Лист2").Cells(7, 17) = SAL
' Расчет относительной ошибки оценки извлекаемых запасов в зависимости от числа лет ' разработки DeltaSAL = Abs((SAL - S_) / S_) Worksheets("Лист2").Cells(34 + n, 1) = SAL Worksheets("Лист2").Cells(34 + n, 2) = DeltaSAL Worksheets("Лист2").Cells(34 + n, 3) = n Worksheets("Лист2").Cells(35, 4) = kr Worksheets("Лист2").Cells(35, 5) = S_
' Расчет прогнозных на L лет значений добычи нефти YL(n + L) = AL(1) * Exp(-AL(2) * (n + L)) * Exp(AL(3) * Log(n + L)) Worksheets("Лист2").Cells(20 + n, 14) = YL(n + L)
' Расчет относительной ошибки прогнозных значений добычи нефти ' DeltaYL1 = ALy(1) * Exp(-ALy(2) * (n + L)) * Exp(ALy(3) * Log(n + L))
DeltaYL1 = y(n) If (DeltaYL1 <> 0) Then DeltaYL = Abs((YL(n) - DeltaYL1) / DeltaYL1) End If Worksheets("Лист2").Cells(20 + n, 15) = DeltaYL1 Worksheets("Лист2").Cells(20 + n, 16) = DeltaYL Worksheets("Лист2").Cells(20 + n, 17) = n Worksheets("Лист2").Cells(21, 18) = L 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 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
End Sub
Function Func1(ByVal t As Single, A() As Single) As Double Func1 = A(1) * ((1 - Exp((-A(2)) * t)) ^ A(3)) End Function
' Процедура определения параметров интгрирвоанной системы моделей добычи нефти методом Гаусса-=Ньютона
Public Sub GN(B1() As Single, C() As Single, Y1() As Single, A() As Single, B() _ As Single, DeltaAL() As Single, m As Integer, ByVal n As Integer, y() As Single, F() As Single, AL() As Single, _ ALh() As Single, kr As Single, MS_ As Single, ByVal t As Integer)
ReDim d(n, m) As Single 'Матрица частных производных(ЧП) в модели объекта ReDim DT(m, n) As Single 'Транспонированная матрица ЧП ReDim DTD(m, m) As Single 'Матрица планирования в методе наименьших квадратов(МНК) ReDim Ds(m, 1) As Single 'Вектор частных производных в модели априорных сведений ReDim DTs(1, m) As Single 'Транспонированная матрица ЧП в модели априорных сведений ReDim DTsDs(m, m) As Single 'Матрица частных производных модели априорной информации ReDim E(n, 1) As Single ' Вектор строка невязок Y-F ReDim DTE(m, 1) As Single ' Вектор правых частей СЛУ в МНК ReDim DTsEs(m, 1) As Single Dim Es As Single ' Es = S_ - SAL Dim QLh(15) As Single ' Функционал качества(для проверки условия сходимости метода ГН) Dim QL(15) As Single Dim i1 As Integer ' Число дробления параметра h Dim i2 As Integer ' Переменная номера шага в методе Гаусса-Нъютона Dim h As Single Dim ALt(3) As Single ReDim FT(n) As Single
h = 1 i2 = 0 Do ' Условный оператор для уточнения параметроа AL(i2+1)=AL(i2)+h*DeltaAL(i2)
' 1. Формирование матрицы частных производных модели добычи нетфи - D(n,m) ' модели объекта f(.,.) For i = 1 To n d(i, 1) = ((1 - Exp((-AL(2)) * i)) ^ AL(3)) d(i, 2) = (i * AL(1) * ((1 - Exp((-AL(2)) * i)) ^ AL(3))) d(i, 3) = (AL(3) * AL(1) * ((1 - Exp((-AL(2)) * i)) ^ (AL(3) - 1))) Next i
' 2. Транспонирование матрицы частных производных - D(n,m)
For i = 1 To m For j = 1 To n DT(i, j) = d(j, i) Next j Next i
' 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) Worksheets("Лист2").Cells(i + 1, 7 + j) = DTD(i, j) Next k Next j Next i
' 4. Формирование вектора Ds модели извлекаемых запасов For i = 1 To m DTs(1, i) = 0 Next i For i = 1 To t DTs(1, 1) = DTs(1, 1) + ((1 - Exp((-AL(2)) * i)) ^ AL(3)) DTs(1, 2) = DTs(1, 2) + (i * AL(1) * ((1 - Exp((-AL(2)) * i)) ^ AL(3))) DTs(1, 3) = DTs(1, 3) + (AL(3) * AL(1) * ((1 - Exp((-AL(2)) * i)) ^ (AL(3) - 1))) Next i
' 5. Транспонирование вектора Ds
For i = 1 To m Ds(i, 1) = DTs(1, i) Worksheets("Лист2").Cells(i + 1, 12) = Ds(i, 1) Next i
' 6. Расчет матрицы DTsDs For i = 1 To m For j = 1 To m DTsDs(i, j) = 0 Next j Next i
For j = 1 To m DTsDs(j, j) = DTsDs(j, j) + Ds(j, 1) * DTs(1, j) Worksheets("Лист2").Cells(j + 1, 13) = DTsDs(j, j) Next j
' 7. Формирование матрицы A=DTD+ krDTsDs
For i = 1 To m For j = 1 To m A(i, j) = DTD(i, j) Next j Next i
For j = 1 To m A(j, j) = A(j, j) + kr * DTsDs(j, j) Next j
For i = 1 To m For j = 1 To m Worksheets("Лист2").Cells(i + 1, 15 + j) = A(i, j) Next j Next i
'8. Формирование вектора столбца Е For i = 1 To n F(i) = Func1(i, AL) E(i, 1) = y(i) - F(i) Next i
'9. Формирование вектора -столбца DTE
For i = 1 To m For k = 1 To n DTE(i, 1) = DT(i, k) * E(k, 1) Worksheets("Лист2").Cells(6 + i, 8) = DTE(i, 1) Next k Next i
'10. Вычисление невязки Es
S = 0 S = S + Func1(t, AL) Es = MS_ - S
'11. Вычисление вектора столбца DTsEs
For i = 1 To m DTsEs(i, 1) = DTs(1, i) * Es Worksheets("Лист2").Cells(6 + i, 12) = DTsEs(i, 1) Next i
'12. Вычисление вектора правых частей ИСМ( B = DTE+kDTsEs) For i = 1 To m B(i) = kr * DTsEs(i, 1) + DTE(i, 1) Worksheets("Лист2").Cells(6 + i, 14) = B(i) Next i
' Решение системы линейных уравнений
Call ELS_hol_1(B1, C, Y1, A, B, DeltaAL, m)
'Вывод решения DeltaAL
For i = 1 To m Worksheets("Лист2").Cells(7 + i2, i + 3) = DeltaAL(i) Next i
'Расчет функционала качества 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 Worksheets("Лист2").Cells(20 + i1 + i2 + i2 * 10, 9) = QLh(i1)
' Определение вектора параметров AL(L)=AL(l-1)+h*Delta AL(L-1) For i = 1 To m ALh(i) = AL(i) Next i
' Проверка условия сходимости метода Гаусса-Ньютона ALt(1) = 178 ALt(2) = 0.21 ALt(3) = 0.65
Do For i = 1 To m AL(i) = ALh(i) + h * DeltaAL(i) Worksheets("Лист2").Cells(20 + i1 + i2 + 10 * i2, i + 5) = AL(i) Worksheets("Лист2").Cells(15, i + 7) = AL(i) Worksheets("Лист1").Cells(13, i + 12) = AL(i) Next i
'Расчет функционала качества QL1 на шаге L
QLh(i1 + 1) = 0
For i = 1 To n F(i) = Func1(i, AL) FT(i) = Func1(i, ALt) QLh(i1 + 1) = QLh(i1 + 1) + (y(i) - F(i)) * (y(i) - F(i)) Next i Worksheets("Лист2").Cells(20 + i1 + i2 + 10 * i2, 10) = QLh(i1 + 1) Worksheets("Лист2").Cells(20 + i1 + i2 + 10 * i2, 11) = i1
' вывод вектора F на шаге i1 For i = 1 To n Worksheets("Лист2").Cells(i + 1, 1) = FT(i) Next i Worksheets("Лист2").Cells(20 + i2 + i1 + 10 * i2, 12) = i2 ' Дробление шага i1 = i1 + 1 h = h / 2 ' Проверка условия сходимости процедуры Гаусса-Нъютона
Loop While (QLh(i1 - 1) >= QLh(1)) And (i1 <= 10)
i2 = i2 + 1 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)
'Окончание процедуры идентификации End Sub
'Вычисление функционала качества Function ArgMin(AL() As Single, ByVal n As Integer, y() As Single, kr As Single, t As Integer) As Single Dim i As Integer Dim S As Single Dim temp As Single Dim SAL As Single S = 0 For i = 1 To n temp = (y(i) - Func1(i, AL)) / y(i) S = S + temp * temp Next i
ArgMin = S
End Function |