Demonstration of Large Scale Model Class LSmodel
This script demonstrates the use of LSmodel class for modeling and model management of Large Scale (LS) systems. The first part shows definition of LSmodel for multiple boilers connected to a single header. The second part is demonstration of LSmodel methods on randomly generated model with 15 submodels.
Contents
- Submodels Definition
- Construction of LSmodel
- Plot model structure
- Random LS model with 15 sub-models
- Structure preserving model order reduction
- "On-line" Model changes
- Submodels grouping
- Time domain models merging
- Overloaded functions
- Impulse response
- Step response
- Bode frequency response
- Nyquist frequency response
- Pole-zero map
- Pole-zero map for I/O pairs
Submodels Definition
This section defines models of individual boilers and header. These models will be used to create LS model in the next section. The structure of the model is in following figure.
% --- Boiler parameters ------------------------------ n_boilers = 5; % number of boilers (1-5) T = [50 100 150 200 250]; % fuel -> steam time constants [s] Ks = 10; % amount of steam [t/hrs] from unit of fuel [t/hrs] V = [300 500 700 800 600]; % boiler volumes [m3] (1/V ~ pressure integration constant) K = [100 130 150 130 140]; % boiler->header pipe "conductivity", =flow/(pressure difference) [t/hrs]/[bar] Kh = 1; % "conductivity" to turbine (higher number -> higher stabilization Vh = 2e3; % header volume [m3] % --- Boiler models ---------------------------------- % Inputs: % FF ... fuel flow - same value flows to all boilers [t/hrs] % ph ... header pressure [MPa] % ---------------------------------------------------- % States: % GS ... generated steam [t/hrs] % pb ... boiler pressure [MPa] % ---------------------------------------------------- % Outputs: % SF ... overall steam flow from all boilers [t/hrs] % pb ... boiler pressure [MPa] % ---------------------------------------------------- for i=1:n_boilers, A = [ -1/T(i) , 0 1/V(i) , -K(i)/V(i) ]; B = [ Ks/T(i) , 0 0 , K(i)/V(i) ]; C = [ 0 , K(i) 0 , 1 ]; D = [ 0 , -K(i) 0 , 0 ]; M{i} = ss(A,B,C,D); % models of individual boilers M{i}.InputName = { 'FF' , 'ph' }; M{i}.StateName = { ['GS' num2str(i)] , ['pb' num2str(i)] }; M{i}.OutputName = { ['SF' num2str(i)] , ['pb' num2str(i)] }; M{i}.Name = ['B' num2str(i) ]; end % --- Header model ----------------------------------- % Inputs: % SF ... total steam flow from boilers [t/hrs] % SD ... steam demand from (from turbine) [t/hrs] % ---------------------------------------------------- % Outputs: % ph ... header pressure [MPa] % ---------------------------------------------------- A = -Kh/Vh; B = [1/Vh -1/Vh]; C = 1; D = 0; M{n_boilers+1} = ss(A,B,C,D); M{n_boilers+1}.InputName = { 'SF' , 'SD' }; M{n_boilers+1}.StateName = { 'ph' }; M{n_boilers+1}.OutputName = { 'ph' }; M{n_boilers+1}.Name = 'H'; % --- Steam flow summing block ----------------------- % sums steam flows from boilers to header sum1 = sumblk('SF','SF1','SF2','SF3','SF4','SF5'); sum1.Name = 'SFsum';
Construction of LSmodel
LS model is contructed by specifying submodels, summators and a list of external inputs and outputs. Internal connections are automatically created by matching input/output names
mod = LSmodel(M,sum1,{'FF','SD'},{'ph','SF'})
Large scale model (total order=11, 6 subsystem(s), 1 summator(s)): --- Subsystems ----------- M1: 2 inputs, 2 outputs, order=2, name="B1" M2: 2 inputs, 2 outputs, order=2, name="B2" M3: 2 inputs, 2 outputs, order=2, name="B3" M4: 2 inputs, 2 outputs, order=2, name="B4" M5: 2 inputs, 2 outputs, order=2, name="B5" M6: 2 inputs, 1 outputs, order=1, name="H" --- Summators ----------------- Sum1: 5 inputs, name = "SFsum" --- External Inputs ----------- IN01: FF X IN02: SD X --- External Outputs ---------- OUT01: ph X OUT02: SF X
Constructor function is compatible with standard CONNECT function, where the same model would be defined as (connect does not accept cell array of models)
mod_connect = connect(M{1},M{2},M{3},M{4},M{5},M{6},sum1,{'FF','SD'},{'ph','SF'});
It is also possible to construct the model by using numeric indexing of inputs and outputs. However, name based signal referencing will be prefered as numeric indexing may be confusing for LS systems.
inputs = { [1 3 5] , 8 }; % extended notation of standard connect
Q = [ 2 7 0 0
4 7 0 0
6 7 0 0
7 1 3 5 ];
outputs = { 7 , [1 3 5] };
inputs = { [1 3 5 7 9], 12 };
Q = [2 11 0 0 0 0
4 11 0 0 0 0
6 11 0 0 0 0
8 11 0 0 0 0
10 11 0 0 0 0
11 1 3 5 7 9 ];
outputs = { 11 , [1 3 5 7 9] };
mod2 = LSmodel(M{1},M{2},M{3},M{4},M{5},M{6},Q,inputs,outputs);
Plot model structure
Plots structure of large scale model as a graph, where vertexes are subsystems and edges indicate interaction between subsystems.
plot(mod);
title('\bfModel Structure');
Computing vertex positions... Done.

Random LS model with 15 sub-models
Previous model is rather small. Larger model will be randomly generated by using command RLS parameterized by number of sub-models, maximum sub-model order, number of external outputs and number of external inputs.
seed=963;rand('seed',seed);randn('seed',seed); % to get stable CL model mod = rls(15,3,3,3)
Large scale model (total order=35, 15 subsystem(s), 8 summator(s)): --- Subsystems ----------- M1: 1 inputs, 1 outputs, order=3, name="subs1" M2: 1 inputs, 1 outputs, order=3, name="subs2" M3: 1 inputs, 2 outputs, order=3, name="subs3" M4: 2 inputs, 2 outputs, order=1, name="subs4" M5: 2 inputs, 1 outputs, order=1, name="subs5" M6: 1 inputs, 1 outputs, order=2, name="subs6" M7: 1 inputs, 1 outputs, order=2, name="subs7" M8: 2 inputs, 2 outputs, order=2, name="subs8" M9: 2 inputs, 2 outputs, order=2, name="subs9" M10: 1 inputs, 1 outputs, order=3, name="subs10" M11: 2 inputs, 2 outputs, order=1, name="subs11" M12: 1 inputs, 2 outputs, order=3, name="subs12" M13: 1 inputs, 2 outputs, order=3, name="subs13" M14: 2 inputs, 2 outputs, order=3, name="subs14" M15: 1 inputs, 1 outputs, order=3, name="subs15" --- Summators ----------------- Sum1: 3 inputs, name = "S1" Sum2: 4 inputs, name = "S2" Sum3: 3 inputs, name = "S3" Sum4: 5 inputs, name = "S4" Sum5: 3 inputs, name = "S5" Sum6: 3 inputs, name = "S6" Sum7: 5 inputs, name = "S7" Sum8: 5 inputs, name = "S8" --- External Inputs ----------- IN01: s24 X IN02: s25 X IN03: s26 X --- External Outputs ---------- OUT01: s7 X OUT02: s8 X OUT03: s13 X
Structure preserving model order reduction
Standard model order reduction destroys model internal structure. Structure preserving model order reduction reduces order of individual sub-models to achieve minimum mismatch with original model from external input/output point of view. Structured singular numbers are plotted to show which sub-models are reduced in order to achieve target global model order.
figure(2);clf; mod_r20=mod.struct_red(20); % model reduced to 20th order mod_r3=mod.struct_red(3); % model reduced to 3rd order % Compare step responses of original and reduced order models figure(1); clf;step(mod,30);hold all; step(mod_r20,'r--',30); step(mod_r3,'g--',30); legend(['Original n=' num2str(mod.n)],['Reduced n=' num2str(mod_r20.n)],['Reduced n=' num2str(mod_r3.n)]);


"On-line" Model changes
Removing submodel(s)
Remove subsystems, which have zero order after structured reduction (just for demonstration of large scale model manipulations).
ind = find(mod_r20.orders==0); removed_models = mod.M(ind); modX = mod.rem_mod(ind); figure(10); clf; subfigure(10,3,3,[1 5]); subplot(121);plot(mod); title('\bfOriginal Model'); subplot(122);plot(modX); title('\bfReduced Model without zero order subsystems');
Computing vertex positions... Done. Computing vertex positions... Done.

Adding submodel(s)
modX = modX.add_mod( removed_models );
cla; plot(modX);
title('\bfRestored Model (different numbering)');
Computing vertex positions... Done.

Removing external input
modX = modX.rem_ext_inp('s24') cla; plot(modX); title('\bfExternal Input Removed');
Large scale model (total order=35, 15 subsystem(s), 8 summator(s)): --- Subsystems ----------- M1: 1 inputs, 1 outputs, order=3, name="subs1" M2: 2 inputs, 2 outputs, order=1, name="subs4" M3: 2 inputs, 1 outputs, order=1, name="subs5" M4: 1 inputs, 1 outputs, order=2, name="subs6" M5: 2 inputs, 2 outputs, order=2, name="subs9" M6: 1 inputs, 1 outputs, order=3, name="subs10" M7: 2 inputs, 2 outputs, order=1, name="subs11" M8: 1 inputs, 2 outputs, order=3, name="subs12" M9: 2 inputs, 2 outputs, order=3, name="subs14" M10: 1 inputs, 1 outputs, order=3, name="subs15" M11: 1 inputs, 1 outputs, order=3, name="subs2" M12: 1 inputs, 2 outputs, order=3, name="subs3" M13: 1 inputs, 1 outputs, order=2, name="subs7" M14: 2 inputs, 2 outputs, order=2, name="subs8" M15: 1 inputs, 2 outputs, order=3, name="subs13" --- Summators ----------------- Sum1: 3 inputs, name = "S1" Sum2: 4 inputs, name = "S2" Sum3: 3 inputs, name = "S3" Sum4: 5 inputs, name = "S4" Sum5: 3 inputs, name = "S5" Sum6: 3 inputs, name = "S6" Sum7: 5 inputs, name = "S7" Sum8: 5 inputs, name = "S8" --- External Inputs ----------- IN01: s25 X IN02: s26 X --- External Outputs ---------- OUT01: s7 X OUT02: s8 X OUT03: s13 X Computing vertex positions... Done.

Adding external input
modX = modX.add_ext_inp('s24') cla; plot(modX); title('\bfExternal Input Added');
Large scale model (total order=35, 15 subsystem(s), 8 summator(s)): --- Subsystems ----------- M1: 1 inputs, 1 outputs, order=3, name="subs1" M2: 2 inputs, 2 outputs, order=1, name="subs4" M3: 2 inputs, 1 outputs, order=1, name="subs5" M4: 1 inputs, 1 outputs, order=2, name="subs6" M5: 2 inputs, 2 outputs, order=2, name="subs9" M6: 1 inputs, 1 outputs, order=3, name="subs10" M7: 2 inputs, 2 outputs, order=1, name="subs11" M8: 1 inputs, 2 outputs, order=3, name="subs12" M9: 2 inputs, 2 outputs, order=3, name="subs14" M10: 1 inputs, 1 outputs, order=3, name="subs15" M11: 1 inputs, 1 outputs, order=3, name="subs2" M12: 1 inputs, 2 outputs, order=3, name="subs3" M13: 1 inputs, 1 outputs, order=2, name="subs7" M14: 2 inputs, 2 outputs, order=2, name="subs8" M15: 1 inputs, 2 outputs, order=3, name="subs13" --- Summators ----------------- Sum1: 3 inputs, name = "S1" Sum2: 4 inputs, name = "S2" Sum3: 3 inputs, name = "S3" Sum4: 5 inputs, name = "S4" Sum5: 3 inputs, name = "S5" Sum6: 3 inputs, name = "S6" Sum7: 5 inputs, name = "S7" Sum8: 5 inputs, name = "S8" --- External Inputs ----------- IN01: s25 X IN02: s26 X IN03: s24 X --- External Outputs ---------- OUT01: s7 X OUT02: s8 X OUT03: s13 X Computing vertex positions... Done.

Removing external output
modX = modX.rem_ext_out('s8') cla; plot(modX); title('\bfExternal Output Removed');
Large scale model (total order=35, 15 subsystem(s), 8 summator(s)): --- Subsystems ----------- M1: 1 inputs, 1 outputs, order=3, name="subs1" M2: 2 inputs, 2 outputs, order=1, name="subs4" M3: 2 inputs, 1 outputs, order=1, name="subs5" M4: 1 inputs, 1 outputs, order=2, name="subs6" M5: 2 inputs, 2 outputs, order=2, name="subs9" M6: 1 inputs, 1 outputs, order=3, name="subs10" M7: 2 inputs, 2 outputs, order=1, name="subs11" M8: 1 inputs, 2 outputs, order=3, name="subs12" M9: 2 inputs, 2 outputs, order=3, name="subs14" M10: 1 inputs, 1 outputs, order=3, name="subs15" M11: 1 inputs, 1 outputs, order=3, name="subs2" M12: 1 inputs, 2 outputs, order=3, name="subs3" M13: 1 inputs, 1 outputs, order=2, name="subs7" M14: 2 inputs, 2 outputs, order=2, name="subs8" M15: 1 inputs, 2 outputs, order=3, name="subs13" --- Summators ----------------- Sum1: 3 inputs, name = "S1" Sum2: 4 inputs, name = "S2" Sum3: 3 inputs, name = "S3" Sum4: 5 inputs, name = "S4" Sum5: 3 inputs, name = "S5" Sum6: 3 inputs, name = "S6" Sum7: 5 inputs, name = "S7" Sum8: 5 inputs, name = "S8" --- External Inputs ----------- IN01: s25 X IN02: s26 X IN03: s24 X --- External Outputs ---------- OUT01: s7 X OUT02: s13 X Computing vertex positions... Done.

Adding external output
modX = modX.add_ext_out('s8') cla; plot(modX); title('\bfExternal Output Added');
Large scale model (total order=35, 15 subsystem(s), 8 summator(s)): --- Subsystems ----------- M1: 1 inputs, 1 outputs, order=3, name="subs1" M2: 2 inputs, 2 outputs, order=1, name="subs4" M3: 2 inputs, 1 outputs, order=1, name="subs5" M4: 1 inputs, 1 outputs, order=2, name="subs6" M5: 2 inputs, 2 outputs, order=2, name="subs9" M6: 1 inputs, 1 outputs, order=3, name="subs10" M7: 2 inputs, 2 outputs, order=1, name="subs11" M8: 1 inputs, 2 outputs, order=3, name="subs12" M9: 2 inputs, 2 outputs, order=3, name="subs14" M10: 1 inputs, 1 outputs, order=3, name="subs15" M11: 1 inputs, 1 outputs, order=3, name="subs2" M12: 1 inputs, 2 outputs, order=3, name="subs3" M13: 1 inputs, 1 outputs, order=2, name="subs7" M14: 2 inputs, 2 outputs, order=2, name="subs8" M15: 1 inputs, 2 outputs, order=3, name="subs13" --- Summators ----------------- Sum1: 3 inputs, name = "S1" Sum2: 4 inputs, name = "S2" Sum3: 3 inputs, name = "S3" Sum4: 5 inputs, name = "S4" Sum5: 3 inputs, name = "S5" Sum6: 3 inputs, name = "S6" Sum7: 5 inputs, name = "S7" Sum8: 5 inputs, name = "S8" --- External Inputs ----------- IN01: s25 X IN02: s26 X IN03: s24 X --- External Outputs ---------- OUT01: s7 X OUT02: s13 X OUT03: s8 X Computing vertex positions... Done.

Extraction of model part - eliminates submodels which are not controlable and observable for selected inputs and outputs.
ex_mod = mod.select('s24','s7'); clf;step(ex_mod); title('\bfStep Response for Selected Model Part');

Submodels grouping
Each sub-model can be assigned to a group. Each group has its own number and sub-models are assigned to given groups in property GROUPING. Following command combines sub-models to groups according to their order.
mod.grouping = mod.orders; clf; plot(mod);
Computing vertex positions... Done.

Time domain models merging
Merging of multiple models describing the same system. The merging is done in time domain and allows to combine information from uncertain ARX models of different orders (idarx or idpoly)
[id_mod,sys] = demo_models; % 3 models of the same system with different % different orders identified from different data for i=1:3, LS{i} = LSmodel(id_mod{i},'u1','y1'); end LS12 = LS{1}.merge('M1',id_mod{2},4,4); LS123 = LS12.merge('M1',id_mod{3},4,4); % --- Plot the result --------------------------- sc = 3; % variance scale figure(10);clf; w = logspace(-1,log10(pi/id_mod{1}.Ts),300); X = [w w(end:-1:1)]; axx = [w(1) w(end) 0 2.5]; for i=1:3, [meanG{i} ,varG{i} ] = LS{i}.freq_uncert('M1',w); end [meanG123,varG123] = LS123.freq_uncert('M1',w); line_c = { 'b','g','m' }; fill_c = { [.8 .8 1],[.8 1 .8],[1 .8 .8]}; for i=1:3, subplot(2,2,i); h(i) = semilogx(w,meanG{i},line_c{i},'LineWidth',3); hold on; Y = meanG{i}+sc*sqrt(varG{i}); Y = [meanG{i}-sc*sqrt(varG{i}) Y(end:-1:1)]; hf(i) = fill(X,Y,fill_c{i},'EdgeAlpha',0,'FaceAlpha',0.5,'LineStyle','none'); semilogx(w,[meanG{i}-sc*sqrt(varG{i}) ; meanG{i}+sc*sqrt(varG{i})]',[line_c{i} '--'],'LineWidth',1); switch i, case 1, title('\bf Model 1 (1st resonance peak)'); case 2, title('\bf Model 2 (2nd resonance peak)'); case 3, title('\bf Model 3 (low frequencies)'); end end subplot(2,2,4); h(4) = semilogx(w,meanG123,'r','LineWidth',3); hold on; Y = meanG123+sc*sqrt(varG123); Y = [meanG123-sc*sqrt(varG123) Y(end:-1:1)]; hf(4) = fill(X,Y,[1 .8 .8],'EdgeAlpha',0,'FaceAlpha',0.5,'LineStyle','none'); semilogx(w,[meanG123-sc*sqrt(varG123) ; meanG123+sc*sqrt(varG123)]','r--','LineWidth',1); title('\bf Merged Models 1, 2, 3'); [mag,phase] = bode(sys,w); for i=1:4, subplot(2,2,i); horig(i) = semilogx(w,squeeze(mag),'k--','LineWidth',3); axis(axx); grid on; xlabel('Frequency (rad/sec)'); ylabel('Magnitude'); legend([h(i) hf(i) horig(i)] ,'mean',['+/- ',num2str(sc),'\sigma'],'orig.'); end
Current iteration error = 30568.7304 Current iteration error = 30568.7304 Current iteration error = 30568.7304 Current iteration error = 370494.737 Current iteration error = 410827.2347 Current iteration error = 68383.1445 Current iteration error = 30568.7304 Current iteration error = 30568.7304 Current iteration error = 30568.7304 Current iteration error = 203078.9615 Maximum number of iterations exceeded. Precission matrix quadratic error (normalized)= 0.02566 Current iteration error = 75719.3471 Current iteration error = 914945.2182 Current iteration error = 870103.539 Current iteration error = 54613.5888 Current iteration error = 155621.7346 Current iteration error = 92255.2874 Current iteration error = 137884.8909 Current iteration error = 863230.4205 Current iteration error = 66342.421 Current iteration error = 947935.3184 Maximum number of iterations exceeded. Precission matrix quadratic error (normalized)= 0.068154 Current iteration error = 4469951.5209 Current iteration error = 2424856.2825 Current iteration error = 2315073.8367 Current iteration error = 2853276.1104 Current iteration error = 2511910.3461 Current iteration error = 2487229.5931 Current iteration error = 2455166.663 Current iteration error = 2764858.2822 Current iteration error = 2181018.6175 Current iteration error = 2194488.5943 Maximum number of iterations exceeded. Precission matrix quadratic error (normalized)= 1.1806 ********** Warning: Fictive data quality may be low. ********** Current iteration error = 2.4564e-020 Precission matrix quadratic error (normalized)= 1.4569e-025

Overloaded functions
Many standard functions are overloaded for LSmodel class.
dcgain(mod)
Warning: Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN. ans = NaN NaN NaN NaN NaN NaN NaN NaN NaN
pole(mod)
ans = -6.6715 -4.4112 + 1.3459i -4.4112 - 1.3459i -0.2786 + 1.6244i -0.2786 - 1.6244i -2.2857 -2.1010 + 0.2367i -2.1010 - 0.2367i -0.1335 + 0.1818i -0.1335 - 0.1818i -0.3904 -1.8142 -0.7314 -0.7882 -0.9910 -1.6759 -1.1631 -1.1784 -1.5414 -1.3276 -1.3349 -1.4047 -0.5040 + 2.7656i -0.5040 - 2.7656i -0.3553 -6.6356 -0.3594 -3.4410 -3.8862 -3.7298 -1.2730 -3.6092 -9.9454 -2.3938 -0.7164
zero(mod)
ans = -7.3180 -6.2397 -0.5040 + 2.7655i -0.5040 - 2.7655i 1.1169 -0.0229 + 1.2125i -0.0229 - 1.2125i -2.8319 -0.0010 -0.8764 + 0.6952i -0.8764 - 0.6952i -0.3551 -0.4630 -0.7853 -2.0245 -1.9684 -1.9588 -1.1695 -1.7014 -1.5418 -1.4130 -1.4413 -6.6356 -0.3594 -3.4410 -3.8862 -3.7298 -1.2730 -3.6092 -9.9454 -2.3938 -0.7164
Impulse response
figure(10);clf; impulse(mod);

Step response
step(mod);

Bode frequency response
bode(mod);

Nyquist frequency response
nyquist(mod);

Pole-zero map
pzmap(mod);

Pole-zero map for I/O pairs
iopzmap(mod);
