-
Notifications
You must be signed in to change notification settings - Fork 0
/
cgClosedLoopOptimization.asv
184 lines (146 loc) · 6.23 KB
/
cgClosedLoopOptimization.asv
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
clear; clc; close all;
mrstModule add ad-core ad-blackoil deckformat ...
agglom upscaling coarsegrid book ...
mrst-gui ad-props incomp optimization...
network-models test-suite linearsolvers co2lab
%% Reload Base Model
%Get Fine Model
load("cgOri_problem.mat")
load("fineRef_data.mat")
disp(cgOri_problem)
%clearPackedSimulatorOutput(fineRef_problem, 'prompt', true) % Remove results
[wellSolCG, statesCG] = getPackedSimulatorOutput(cgOri_problem);
%Retrieve Data
cgModel = cgOri_problem.SimulatorSetup.model;
cgSchedule = cgOri_problem.SimulatorSetup.schedule;
W_cg = cgOri_problem.SimulatorSetup.schedule.control.W;
ts_cg = cgOri_problem.SimulatorSetup.schedule.step.val;
G_cg = cgOri_problem.SimulatorSetup.model.G;
%upscale initial schedule
cState0 = upscaleState(cgModel, model, initState);
%case Settings
baseCaseName = 'optim_test_NPVwithContainment';
fine_baseName= 'cgModelOri_';
reset_flag = true; %true to delete prev problem run
plotToolbar(G_cg,statesCG)
plotWell(G_cg,W_cg)
plotWellSols(wellSolCG)
%% Define Containment Region
% Penalty Area:
K_target = 7;
mon_idx = false(G_cg.parent.cartDims);
%mon_idx(32:29,18:89,K_target) = true; %I1-P1
mon_idx(39:47,2:60,K_target) = true;
%mon_idx(22:55,65:87,K_target)= true;
cid_mon = find(mon_idx(G_cg.parent.cells.indexMap));
cid_mon_CG = cgModel.G.partition(cid_mon);
plotGrid(cgModel.G, 'facecolor', 'none', 'edgealpha', 0.1);
plotGrid(cgModel.G,cid_mon_CG,'facecolor','red','edgealpha',0.1);
plotWell(cgModel.G,W_cg,'color','blue','color2','red')
%% Fine Scale Full Optimization
%% Model Optimization: NPV CO2
%Constraints
bhp1 = W_cg(3).val;
li = [0.01 10]* mega * 1e3 / year / cgModel.fluid.rhoGS; % Injector limits (MT/yr)
lp = [0.1*bhp1 bhp1]; % Producer limits
scaling=[];
scaling.boxLims = [li;li;lp;lp]; % control scaling
scaling.obj = 3.2e7; % objective scaling
% Get initial scaled controls
% u_base = schedule2control(schedule_pred, scaling); %add pred schedule restart
% Define objective function
d = 0.05; % yearly discount factor
ro = 150; % co2 produced handling cost ($/ton)
rwp = 6; % water production handling costs ($/stb)
rwi = 87; % carbon price ($/ton)
npvopts = {'CarbonProductionCost', ro , ...
'WaterProductionCost', rwp , ...
'CarbonPrice', rwi , ...
'DiscountFactor', d};
npvFn = @(model, states, schedule, varargin)NPVCO2(model, states, schedule,varargin{:}, npvopts{:});
%% Forward Optimization
% Optimize next control steps
timePerOpt = 50;
ControlPerYear = 1;
%% ---set forecast schedule [1 control per year with 12 ts(mon)/ctrl]
ts = transpose(repmat(1/1, 1, timePerOpt*ControlPerYear)*year); %repmat(1/12, 12, timePerOpt*ControlPerYear)*year
ts = transpose(mat2cell(ts,ones(timePerOpt*ControlPerYear,1)));
cg_sched_base = [];
numCnt = numel(ts);
for i=1:numCnt
cg_sched_base.control(1,i).W = W_cg;
end
cg_sched_base.step.control = rldecode((1:numCnt)', cellfun(@numel, ts));
cg_sched_base.step.val = transpose(horzcat(ts{:}));
u_base = schedule2control(cg_sched_base, scaling); %add pred schedule restart
%% ---Get function handle for objectiveclear evaluation
f = @(u)evalObjective(u, npvFn, cState0, cgModel, cg_sched_base, scaling);
%---Optimize control u
%!! TODO Change to Optimization Problem
[vOpt, u_opt, hisOpt] = unitBoxBFGS(u_base, f,'objChangeTol', 1e-8, ... %solving optimization step and return optimal parameter popt
'gradTol', 1e-5, 'maxIt',9, 'lbfgsStrategy', 'dynamic', ...
'lbfgsNum', 5, 'outputHessian',true, 'logPlot', true,'maximize',true);
%save obj fun val and optimizer history
objValOpt=vOpt;
%historyOpt{iter}=hisOpt;
save("historyOpt.mat","hisOpt")
sched_opt = control2schedule(u_opt, cg_sched_base, scaling);
save("sched_opt.mat","sched_opt");
%% Forward Run using Optimized Schedule
fineOpt_problem = packSimulationProblem(cState0,cgModel,sched_opt,...
baseCaseName,'Name','fine_opt_fullwindow',...
'Description','optimized fine model full window');
if reset_flag==true
clearPackedSimulatorOutput(fineOpt_problem, 'prompt', false) % Remove results
disp('previous results deleted')
end
save("fineOptProblem_.mat","fineOpt_problem") %save problem
simulatePackedProblem(fineOpt_problem);
[wellSolFineOpt,statesFineOpt] = getPackedSimulatorOutput(fineOpt_problem);
%% Forward run of Base Model
cgSchedule.step.val = ones(timePerOpt,1)*year;
cgSchedule.step.val = ones(timePerOpt,1)*year;
%% Plot Well Data
fh_well = plotWellSols({wellSolCG(1:timePerOpt),wellSolFineOpt},...
{cgSchedule.step.val(1:timePerOpt),sched_opt.step.val},...
'datasetnames',{'basecase','optimized'}, ...
'zoom', false, ...
'timescale','years',...
'field', 'qGs');
%% Plot Schedules
figure
plotSchedules(sched_opt, 'singlePlot', true, 'boxConst', [li;li;lp;lp] )
figure
plotSchedules(cg_sched_base, 'singlePlot', true, 'boxConst', [li;li;lp;lp] )
%% Plot States
figure
plotToolbar(G_cg,statesFineOpt);
plotWell(G_cg,W_cg);
figure
plotToolbar(G_cg,statesCG(1:timePerOpt));
plotWell(G_cg,W_cg);
%%
schedFullRef.control.W = cg_sched_base.control.W;
schedFullRef.step.val = cg_sched_base.step.val(1:timePerOpt);
schedFullRef.step.control = cg_sched_base.step.control(1:timePerOpt);
load("sched_opt.mat")
vals = cell2mat(NPVCO2(cgModel, statesCG(1:timePerOpt), schedFullRef, npvopts{:}));
vals_opt = cell2mat(NPVCO2(cgModel, statesFineOpt, sched_opt, npvopts{:}));
dtime = schedFullRef.step.val/year;
time = cumsum(dtime);
% Plot discounted net cashflow $/day:
figure, plot(time, vals./dtime, '--b','LineWidth', 2);
hold on, plot(time, vals_opt./dtime, '-b','LineWidth', 2);
line([0 30], [0 0], 'color', 'r'), set(gca, 'FontSize', 14)
title('Net cash-flow [$]'), legend('Base', 'Optimal')
% Find index of first occuring time > 10 days, where net cashflow becomes
% negative:
inx = find(and(vals<0, time>10), 1, 'first');
% Plot evolution of NPV and indicate peak value:
npv = cumsum(vals);
figure, plot(time, cumsum(vals), '--b', 'LineWidth', 2);
hold on, plot(time, cumsum(vals_opt), '-b', 'LineWidth', 2);
%plot([1 1]*time(inx), [0 npv(inx)], '--k', 'LineWidth', 2)
set(gca, 'FontSize', 14), title('Evolution of NPV [$]'),
legend('Base', 'Optimal', 'Location', 'northwest')
hold off