-
Notifications
You must be signed in to change notification settings - Fork 2
/
mpod_ro.m
49 lines (36 loc) · 1.7 KB
/
mpod_ro.m
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
function [proj,name] = eds_ro(solver,discrete,scenario,config)
%%% project: morgen - Model Order Reduction for Gas and Energy Networks
%%% version: 1.2 (2022-10-07)
%%% authors: C. Himpe (0000-0003-2194-6754), S. Grundel (0000-0002-0209-6566)
%%% license: BSD-2-Clause (opensource.org/licenses/BSD-2-clause)
%%% summary: Structured nonlinear modified POD.
global ODE;
name = 'Modified POD RO';
logger('head',name);
iP = 1:discrete.nP;
iQ = discrete.nP+1:discrete.nP+discrete.nQ;
sysdim = [discrete.nPorts, discrete.nP + discrete.nQ, discrete.nPorts];
timedisc = [config.solver.dt, scenario.tH];
flags = [3,0,0,1,1,0,0,0,0,0,0,0,0];
% Specialize state solver
ODE = @(f,g,t,x0,u,p) solver(setfields(discrete,'C',1,'x0',x0), ...
setfields(scenario,'ut',u,'T0',p(1),'Rs',p(2)), ...
config.solver).y;
% Empirical reachability Gramian
WR = emgr(@() 0,@() 1,sysdim,timedisc,'c',config.samples,flags,config.excitation,[],[],0.01*scenario.us);
% Specialize output solver
ODE = @(f,g,t,x0,u,p) solver(setfields(discrete,'x0',x0), ...
setfields(scenario,'ut',u,'T0',p(1),'Rs',p(2)), ...
config.solver).y;
% Empirical observability Gramian
WO = emgr(@() 0,@(x,u,p,t) discrete.C * x,sysdim,timedisc,'o',config.samples,flags,config.excitation);
% Pressure projector
[pc,sc,~] = svds(WR(iP,iP),config.rom_max);
[po,so,~] = svds(WO(iP,iP),config.rom_max);
% Mass-flux projector
[qc,sc,~] = svds(WR(iQ,iQ),config.rom_max);
[qo,so,~] = svds(WO(iQ,iQ),config.rom_max);
proj = {po,pc; ...
qo,qc};
ODE = [];
end