balancmr

If you do not specify any target order for the reduced model, balancmr displays the Hankel singular values of the model and prompts you to choose a reduced-model order.

For this example, use a random 30th-order state-space model.

rng(1234,'twister');     % fix random seed for example repeatability
G = rss(30,5,4);

G1 = balancmr(G)

Please enter the desired order: (>=0)

sigma(G,G1)

When you have particular target order or orders in mind, you can use balancmr to reduce a high-order model to those orders. For this example, use a random 30th-order state-space model.

rng(1234,'twister');     % fix random seed for example repeatability
G = rss(30,5,4);

Use a scalar input argument to reduce the model to a single order. For example, compute a 20th-order approximation.

[G1,info1] = balancmr(G,20);
sigma(G,G1)

Use a vector to generate several approximations. The following command returns an array of models of even orders from 10 to 18.

[G2,info2] = balancmr(G,[10:2:18]);
sigma(G,G2)

Obtain the lowest-order approximation such that the sum of the Hankel singular values of the truncated states does not exceed a specified value. For this example, use a random 30th-order state-space model.

rng(1234,'twister');     % fix random seed for example repeatability
G = rss(30,5,4);

Compute two approximate models, one for which the error does not exceed 0.1, and a second for which the error does not exceed 0.5. To do so, provide these values in an array. balancmr returns an array of approximate models.

Gr = balancmr(G,'MaxError',[0.1 0.5]);
size(Gr)

2x1 array of state-space models.
Each model has 5 outputs, 4 inputs, and between 27 and 29 states.

Examine the results.

sigma(G,Gr)

Reduce a 4th-order system to a second-order approximation with emphasis on the frequency band 10 rad/s - 100 rad/s. Consider the following system.

sys = tf(1,[1 0.5 1]) + tf(100*[1/10 1],[1 10 1000]); 
bode(sys)

To focus the model-reduction algorithm on the higher-frequency dynamics, specify a function with a bandpass profile.

s = tf('s'); 
w1 = (s+1)/(s/10+1)/(s/60+1)*(s/600+1); 
bodemag(w1)

The plot confirms that the weighting function w1 has the desired profile, peaking between 10 rad/s and 100 rad/s. To perform the reduction, specify the inverse of this profile as the output weight, using the 'Weights' option of balancmr.

weight = {1/w1,1}; 
wrsys = balancmr(sys,2,'Weights',weight);

Compare the result with a second-order model obtained without the weighting.

rsys = balancmr(sys,2);
bode(sys,rsys,wrsys) 
legend('Original','Unweighted','Weighted')

The model obtained with the weighting function provides a better match for the dynamics in the frequency band 10 rad/s - 100 rad/s.

Use a weighting function to control the frequency dependence of the error between the original and reduced models.

For this example, load a 48-state SISO model and reduce it to 6th order.

load('balancmrData.mat','bplant')
bplant6 = balancmr(bplant,6);
bode(bplant,bplant6,bplant-bplant6,logspace(0,2,200))
legend('Original','Reduced','Error')

The error is fairly uniformly distributed in frequency. Create a weighting function that allows for a larger error at frequencies below 100 rad/s. In this case, use a biproper function that has unit gain at low frequency, but drops to -40 dB at higher frequencies.

W = tf(0.01*[1 1.4e2 1e4],[1 14 100]); 
bodemag(W,bplant)

Reduce the model to 6th order again, using this weighting function. Because bplant is a SISO model, you can use the function as either input or output weight.

bplant6W = balancmr(bplant,6,'Weights',{W,1});
bode(bplant,bplant6W,bplant-bplant6W,logspace(0,2,200))
legend('Original','Reduced w/Weight','Error')

The error is now larger at low frequencies, with a correspondingly better match at high frequencies between the original and reduced models.