simplifying higher-凯发k8网页登录

introduction

robust control toolbox offers tools to deal with large models such as:

for control purposes, it is generally enough to have an accurate model near the crossover frequency. for simulation, it is enough to capture the essential dynamics in the frequency range of the excitation signals. this means that it is often possible to find low-order approximations of high-order models. robust control toolbox offers a variety of model-reduction algorithms to best suit your requirements and your model characteristics.

the model reduction process

a model reduction task typically involves the following steps:

we can easily access these algorithms through the top-level interface reduce. the methods employ different measures of "closeness" between the original and reduced models. the choice is application-dependent. let's try each of them to investigate their relative merits.

example: a model for rigid body motion of a building

in this example, we apply the reduction methods to a model of the building of the los angeles university hospital. the model is taken from slicot working note 2002-2, "a collection of benchmark examples for model reduction of linear time invariant dynamical systems," by y. chahlaoui and p.v. dooren. it has eight floors, each with three degrees of freedom - two displacements and one rotation. we represent the input-output relationship for any one of these displacements using a 48-state model, where each state represents a displacement or its rate of change (velocity).

let's load the data for the example:

load buildingdata.mat

examining the plant dynamics

let's begin by analyzing the frequency response of the model:

bode(g)
grid on

figure 1: bode diagram to analyze the frequency response

as observed from the frequency response of the model, the essential dynamics of the system lie in the frequency range of 3 to 50 radians/second. the magnitude drops in both the very low and the high-frequency ranges. our objective is to find a low-order model that preserves the information content in this frequency range to an acceptable level of accuracy.

computing hankel singular values

to understand which states of the model can be safely discarded, look at the hankel singular values of the model:

hsv_add = hankelsv(g);
bar(hsv_add)
title('hankel singular values of the model (g)');
xlabel('number of states')
ylabel('singular values (\sigma_i)')
line([10.5 10.5],[0 1.5e-3],'color','r','linestyle','--','linewidth',1)
text(6,1.6e-3,'10 dominant states.')

figure 2: hankel singular values of the model (g).

the hankel singular value plot suggests that there are four dominant modes in this system. however, the contribution of the remaining modes is still significant. we'll draw the line at 10 states and discard the remaining ones to find a 10th-order reduced model gr that best approximates the original system g.

performing model reduction using an additive error bound

the function reduce is the gateway to all model reduction routines available in the toolbox. we'll use the default, square-root balance truncation ('balancmr') option of reduce as the first step. this method uses an "additive" error bound for reduction, meaning that it tries to keep the absolute approximation error uniformly small across frequencies.

% compute 10th-order reduced model (reduce uses balancmr method by default)
[gr_add,info_add] = reduce(g,10);
% now compare the original model g to the reduced model gr_add
bode(g,'b',gr_add,'r')
grid on
title('comparing original (g) to the reduced model (gr\_add)')
legend('g - 48-state original ','gr\_add - 10-state reduced','location','northeast')

figure 3: comparing original (g) to the reduced model (gr_add)

performing model reduction using a multiplicative error bound

as seen from the bode diagram in figure 3, the reduced model captures the resonances below 30 rad/s quite well, but the match in the low frequency region (<2 rad/s) is poor. also, the reduced model does not fully capture the dynamics in the 30-50 rad/s frequency range. a possible explanation for large errors at low frequencies is the relatively low gain of the model at these frequencies. consequently, even large errors at these frequencies contribute little to the overall error.

to get around this problem, we can try a multiplicative-error method such as bstmr. this algorithm emphasizes relative errors rather than absolute ones. because relative comparisons do not work when gains are close to zero, we need to add a minimum-gain threshold, for example by adding a feedthrough gain d to our original model. assuming we are not concerned about errors at gains below -100 db, we can set the feedthrough to 1e-5.

gg = g;
gg.d = 1e-5;

now, let's look at the singular values for multiplicative (relative) errors (using the 'mult' option of hankelsv)

hsv_mult = hankelsv(gg,'mult');
bar(hsv_mult)
title('multiplicative-error singular values of the model (g)');
xlabel('number of states')
ylabel('singular values (\sigma_i)')

figure 4: multiplicative-error singular values of the model (g)

a 26th-order model looks promising, but for the sake of comparison to the previous result, let's stick to a 10th order reduction.

% use bstmr algorithm option for model reduction 
[gr_mult,info_mult] = reduce(gg,10,'algorithm','bst');
%now compare the original model g to the reduced model gr_mult
bode(g,gr_add,gr_mult,{1e-2,1e4}), grid on
title('comparing original (g) to the reduced models (gr\_add and gr\_mult)')
legend('g - 48-state original ','gr\_add (balancmr)','gr\_mult (bstmr)','location','northeast')

figure 5: comparing original (g) to the reduced models (gr_add and gr_mult)

the fit between the original and the reduced models is much better with the multiplicative-error approach, even at low frequencies. we can confirm this by comparing the step responses:

step(g,gr_add,gr_mult,15) %step response until 15 seconds
legend('g: 48-state original ','gr\_add: 10-state (balancmr)','gr\_mult: 10-state (bstmr)')

figure 6: step responses of the three models

validating the results

all algorithms provide bounds on the approximation error. for additive-error methods like balancmr, the approximation error is measured by the peak (maximum) gain of the error model g-greduced across all frequencies. this peak gain is also known as the h-infinity norm of the error model. the error bound for additive-error algorithms looks like:

where the sum is over all discarded hankel singular values of g (entries 11 through 48 of hsv_add). we can verify that this bound is satisfied by comparing the two sides of this inequality:

norm(g-gr_add,inf) % actual error

ans = 6.0251e-04

% theoretical bound (stored in the "errorbound" field of the "info" 
% struct returned by |reduce|)
info_add.errorbound

ans = 0.0047

for multiplicative-error methods such as bstmr, the approximation error is measured by the peak gain across frequency of the relative error model g\(g-greduced). the error bound looks like

where the sum is over the discarded multiplicative hankel singular values computed by hankelsv(g,'mult'). again we can compare these bounds for the reduced model gr_mult

norm(gg\(gg-gr_mult),inf) % actual error

ans = 0.5949

% theoretical bound 
info_mult.errorbound

ans = 546.1730

plot the relative error for confirmation

bodemag(gg\(gg-gr_mult),{1e-2,1e3})
grid on
text(0.1,-50,'peak gain: -4.6 db (59%) at 17.2 rad/s')
title('relative error between original model (g) and reduced model (gr\_mult)')

figure 7: relative error between original model (g) and reduced model (gr_mult)

from the relative error plot above, there is up to 59% relative error at 17.2 rad/s, which may be more than we're willing to accept.

picking the lowest order compatible with a desired accuracy level

to improve the accuracy of gr_mult, we'll need to increase the order. to achieve at least 5% relative accuracy, what is the lowest order we can get? the function reduce can automatically select the lowest-order model compatible with our desired level of accuracy.

% specify a maximum of 5% approximation error
[gred,info] = reduce(gg,'errortype','mult','maxerror',0.05);
size(gred)

state-space model with 1 outputs, 1 inputs, and 35 states.

the algorithm has picked a 34-state reduced model gred. compare the actual error with the theoretical bound:

norm(gg\(gg-gred),inf)

ans = 0.0068

info.errorbound

ans = 0.0342

look at the relative error magnitude as a function of frequency. higher accuracy has been achieved at the expense of a larger model order (from 10 to 34). note that the actual maximum error is 0.6%, much less than the 5% target. this discrepancy is a result of the function bstmr using the error bound rather than the actual error to select the order.

bodemag(gg\(gg-gred),{1,1e3})
grid on
text(5,-75,'peak gain: -43.3 db (0.6%) at 73.8 rad/s')
title('relative error between original model (g) and reduced model (gred)')

figure 8: relative error between original model (g) and reduced model (gred)

compare the bode responses

bode(g,gred,{1e-2,1e4})
grid on
legend('g - 48-state original','gred - 34-state reduced')

figure 9: bode diagram of the 48-state original model and the 34-state reduced model

finally, compare the step responses of the original and reduced models. they are virtually indistinguishable.

step(g,'b',gred,'r--',15) %step response until 15 seconds
legend('g: 48-state original ','gred: 34-state (bstmr)')
text(5,-4e-4,'maximum relative error <0.05')

figure 10: step response plot of the 48-state original model and the 34-state reduced model