act8pz.m

vibration simulation using ansys and matlab 一书中的程序

	index_h1_coil_sort = fliplr(index_h1_coil_sort)				% max to min indices

	index_h0_piez_sort = fliplr(index_h0_piezo_sort)				% max to min indices

	index_h1_piez_sort = fliplr(index_h1_piezo_sort)				% max to min indices

    index_orig = 1:num_modes_total;

	[index_h0_coil_sort' index_h1_coil_sort' index_h0_piez_sort' index_h1_piez_sort']

%	plot results
 
   	semilogy(index_orig(2:num_modes_total),freqvec(2:num_modes_total),'k-');
   	title(['frequency versus mode number'])
   	xlabel('mode number')
   	ylabel('frequency, hz')
   	grid off
	disp('execution paused to display figure, "enter" to continue'); pause

   	semilogy(index_orig,gain_h0_coil,'k.-',index_orig,gain_h1_coil,'k-')
   	title(['coil input:  dc value of each mode contribution versus mode number'])
   	xlabel('mode number')
   	ylabel('dc value')
	legend('h0 coil input','h1 coil input')
   	grid off
	disp('execution paused to display figure, "enter" to continue'); pause

   	semilogy(index_orig,gain_h0_piezo,'k.-',index_orig,gain_h1_piezo,'k-')
   	title(['piezo input:  dc value of each mode contribution versus mode number'])
   	xlabel('mode number')
   	ylabel('dc value')
	legend('h0 piezo input','h1 piezo input')
   	grid off
	disp('execution paused to display figure, "enter" to continue'); pause

   	loglog(freqvec(2:num_modes_total),gain_h0_coil(2:num_modes_total),'k.-', ...
   		   freqvec(2:num_modes_total),gain_h1_coil(2:num_modes_total),'k-')
   	title(['coil input:  dc value of each mode contribution versus frequency'])
   	xlabel('frequency, hz')
   	ylabel('dc value')
	axis([500 25000 -inf inf])
	legend('h0 coil input','h1 coil input')
   	grid off
	disp('execution paused to display figure, "enter" to continue'); pause

   	loglog(freqvec(2:num_modes_total),gain_h0_piezo(2:num_modes_total),'k.-', ...
   		   freqvec(2:num_modes_total),gain_h1_piezo(2:num_modes_total),'k-')
   	title(['piezo input:  dc value of each mode contribution versus frequency'])
   	xlabel('frequency, hz')
   	ylabel('dc value')
	axis([500 25000 -inf inf])
	legend('h0 piezo input','h1 piezo input')
   	grid off
	disp('execution paused to display figure, "enter" to continue'); pause

   	semilogy(index_orig,gain_h0_coil_sort,'k.-',index_orig,gain_h1_coil_sort,'k-')
   	title(['coil input:  sorted dc value of each mode versus number of modes included'])
   	xlabel('modes included')
   	ylabel('sorted dc value')
	legend('h0 coil input','h1 coil input')
   	grid off
	disp('execution paused to display figure, "enter" to continue'); pause

   	semilogy(index_orig,gain_h0_piezo_sort,'k.-',index_orig,gain_h1_piezo_sort,'k-')
   	title(['piezo input:  sorted dc value of each mode versus number of modes included'])
   	xlabel('modes included')
   	ylabel('sorted dc value')
	legend('h0 piezo input','h1 piezo input')
   	grid off
	disp('execution paused to display figure, "enter" to continue'); pause

   	semilogy(index_orig,gain_h0_coil_sort,'k.-',index_orig,gain_h1_coil_sort,'k.-', ...
   			 index_orig,gain_h0_piezo_sort,'k-',index_orig,gain_h1_piezo_sort,'k-')
   	title(['coil and piezo input:  sorted dc value of each mode versus number of modes included'])
   	xlabel('modes included')
   	ylabel('sorted dc value')
	legend('h0 coil input','h1 coil input','h0 piezo input','h1 piezo input')
   	grid off
	disp('execution paused to display figure, "enter" to continue'); pause

%	create five state space systems with all modes included, differing in the ordering
%	of the modes, the unsorted system will be used for all reductions, letting balreal do all
%	the ordering, the sorted systems will be used to show how the dc gain ordering compares 
%	with the balanced ordering

%		1)  unsorted
%		2)  sorted, head 0, coil input
%		3)  sorted, head 1, coil input
%		4)  sorted, head 0, piezo input
%		5)  sorted, head 1, piezo input

	for  num_model = 1:5

	if  num_model == 1			% unsorted

		xnnew = xn;

		freqnew = freqvec;

	elseif  num_model == 2		% sorted, head 0, coil input

		xnnew = xn(:,index_h0_coil_sort);

		freqnew = freqvec(index_h0_coil_sort);

	elseif  num_model == 3		% sorted, head 1, coil input

		xnnew = xn(:,index_h1_coil_sort);

		freqnew = freqvec(index_h1_coil_sort);

	elseif  num_model == 4		% sorted, head 0, piezo input

		xnnew = xn(:,index_h0_piezo_sort);

		freqnew = freqvec(index_h0_piezo_sort);

	elseif  num_model == 5		% sorted, head 1, piezo input

		xnnew = xn(:,index_h1_piezo_sort);

		freqnew = freqvec(index_h1_piezo_sort);
	
	end

%	define variables for all modes included system matrix, a

    w = freqnew*2*pi;		%  	frequencies in rad/sec

	w2 = w.^2;				

    zw = 2*zeta_unsort.*w;

%	define size of system matrix

	asize = 2*num_modes_total;

   	disp('   ');
   	disp('   ');
   	disp(['size of system matrix a is ',num2str(asize)]);

%  	setup system matrix for all modes included model

   	a = zeros(asize);

   	for  col = 2:2:asize

      	row = col-1;

      	a(row,col) = 1;

   	end

   	for  col = 1:2:asize

      	row = col+1;

      	a(row,col) = -w2((col+1)/2);

   	end

   	for  col = 2:2:asize

      	row = col;

      	a(row,col) = -zw(col/2);

   	end

%	setup input matrix b, state space forcing function in principal coordinates

%	two-input system

%		first input is coil force
%		second input is excitation of both piezo elements with opposite polarity

	f_physical = [f_coil f_piezo];
	
%	f_principal is the matrix of forces in principal coordinates

   	f_principal = xnnew'*f_physical;

%	b is the matrix of forces in principal coordinates, state space form

	b = zeros(2*num_modes_total,2);

	for  cnt = 1:num_modes_total
	
		b(2*cnt,:) = f_principal(cnt,:);
		
	end				

%  	setup cdisp and cvel, padded xn matrices to give the displacement and velocity
%	vectors in physical coordinates
%	cdisp and cvel each have numdof rows and alternating columns consisting of columns
%   of xnnew and zeros to give total columns equal to the number of states

%	all modes included cdisp and cvel

   	for  col = 1:2:2*length(freqnew)

      	for  row = 1:numdof

         	c_disp(row,col) = xnnew(row,ceil(col/2));

         	cvel(row,col) = 0;
   
      	end

   	end

   	for  col = 2:2:2*length(freqnew)

      	for  row = 1:numdof

         	c_disp(row,col) = 0;

	        cvel(row,col) = xnnew(row,col/2);

      	end

   	end

%  	define output

   	d = [0];  %

	if  num_model == 1			% unsorted

		sys = ss(a,b,c_disp(31:32,:),d);

	elseif  num_model == 2		% sorted, head 0, coil input

		sys_h0_coil = ss(a,b,c_disp(31:32,:),d);

	elseif  num_model == 3		% sorted, head 1, coil input

		sys_h1_coil = ss(a,b,c_disp(31:32,:),d);

	elseif  num_model == 4		% sorted, head 0, piezo input

		sys_h0_piezo = ss(a,b,c_disp(31:32,:),d);

	elseif  num_model == 5		% sorted, head 1, piezo input

		sys_h1_piezo = ss(a,b,c_disp(31:32,:),d);

	end

	end   		% end of for loop for creating system matrices

%	partition system matrices into rigid body mode and oscillatory modes, can't use balreal
%	with rigid body mode so will reduce the oscillatory modes and then augment the resulting
%	system with the rigid body mode

%	define oscillatory system, where output 31 is head 1, output 32 is head 0

	[a,b,c_disp,d] = ssdata(sys);

	a_syso = a(3:asize,3:asize);

	b_syso = b(3:asize,:);

	c_syso = c_disp(1:2,3:asize);

	syso = ss(a_syso,b_syso,c_syso,d);

%	define controllability and observability gramians for oscillatory system, syso

	wc = gram(syso,'c');

	wo = gram(syso,'o');

	[row_syso,col_syso] = size(a_syso);

	statevec = 1:row_syso;

%	plot controllability and observability gramians

	meshz(wc);
	view(60,30);
	title(['controllability gramian for oscillatory system'])
	xlabel('state')
	ylabel('state')
	grid on

	disp('execution paused to display figure, "enter" to continue'); pause

	meshz(wo);
	view(60,30);
	title(['observability gramian for oscillatory system'])
	xlabel('state')
	ylabel('state')
	grid on

	disp('execution paused to display figure, "enter" to continue'); pause

%	pull out diagonal elements

	wc_diag = diag(wc);

	wo_diag = diag(wo);

%	plot diagonal terms of controllability and observability gramians

	semilogy(statevec,wc_diag,'k.-')
	title(['controllability gramian diagonal terms'])
	xlabel('states')
	ylabel('diagonal')
	grid off

	disp('execution paused to display figure, "enter" to continue'); pause

	semilogy(statevec,wo_diag,'k.-')
	title(['observability gramian diagonal terms'])
	xlabel('states')
	ylabel('diagonal')
	grid off

	disp('execution paused to display figure, "enter" to continue'); pause

%	position and velocity states plotted separately

	semilogy(statevec(1:2:row_syso),wc_diag(1:2:row_syso),'k.-', ...
		 					statevec(2:2:row_syso),wc_diag(2:2:row_syso),'k-')
	title(['controllability gramian diagonal terms'])
	xlabel('states')
	ylabel('diagonal')
	legend('position states','velocity states',3)
	grid off

	disp('execution paused to display figure, "enter" to continue'); pause

	semilogy(statevec(1:2:row_syso),wo_diag(1:2:row_syso),'k.-', ...
		  	 				statevec(2:2:row_syso),wo_diag(2:2:row_syso),'k-')
	title(['observability gramian diagonal terms'])
	xlabel('states')
	ylabel('diagonal')
	legend('position states','velocity states',3)
	grid off

	disp('execution paused to display figure, "enter" to continue'); pause

%	use balreal to rank oscillatory states and modred to reduce for comparison

	[sysob,g,T,Ti] = balreal(syso);

	[ao_bal,bo_bal,cdispo_bal,do_bal] = ssdata(sysob);

	semilogy(g,'k.-')
	title('diagonal of balanced gramian versus number of states')
	xlabel('state number')
	ylabel('diagonal of balanced gramian')
	grid off

	osc_states_used = input(['enter number of oscillatory states to use, default 20  ...  ']);

	if  isempty(osc_states_used)

		osc_states_used = 20;

	end