Modeling is the art of deriving equations to represent the pertinent.
For proposals, early stage designs, quick turn around analyses a frequency domain model is often superior to a time domain model. Frequency domain models take disturbance PSDs directly, use transfer functions directly, and produce output or residual PSDs directly. The answer is a steady-state response. Often times the controller is shooting for 0 so the steady-state reponse is also the residual error that will be the analysis output or metric for report.
2 Brief Overview of the Math
Frequncy domain modeling is a matter of determining the impulse response of a system to a random process. Typically the theory for this is detailed in texts on stochastic processes and stochastic contol.
|||Freq Domain #1|
- is the one-sided input PSD in
- is the frequency response function of the system and
- is the one-sided output PSD or auto power spectral density function.
The frequency response function, , is related to the impulse response function (transfer function) by
Note some texts will state that this is only valid for random processes which are stationary. Other texts suggest stationary and ergodic while still others state weakly stationary processes. Some texts do not distinguish between strictly stationary and weakly stationary. From practice, the rule thumb is if the PSD of the input process is the same from hour to hour and day to day then the input PSD can be used and the above equation is valid.
3 In Practice
The standard system block diagram shown in Figure 2 represents many systems well.
In order to predict the systems response to an input disturbance we must calculate the system's disturbance rejection transfer function. From the block diagram quick reference the disturbance rejection is
- is the disturbance rejection transfer function
Note that is shorthand for and that the is from to .
The magnitude, at discrete frequencies, of the disturbance rejection transfer function becomes in Freq Domain #1 above.
3.1 MATLAB code
Assume we have a system with a closed loop bandwidth of approximately 1 kHz and little peaking such that
>> wn = 1000*(2*pi); >> z = 0.5; >> k = 1; >> CL = tf([wn^2], [1 2*z*wn wn^2]);
Then the disturbance rejection is
>> DR = 1 - CL;
To get the magnitude of each use the bode command
>> freqVec = logspace(-1, 4, 5000); >> [CLmag, CLphs] = bode(CL, freqVec*(2*pi)); >> [DRmag, DRphs] = bode(DR, freqVec*(2*pi));
Plot the transfer functions to double check them
>> figure >> semilogx(freqVec, db(CLmag(:)), 'b', freqVec, db(DRmag(:)), 'r-.') >> grid on >> xlabel('Frequency (Hz)') >> ylabel('Magnitude (dB)') >> title('Example Closed Loop and Disturbance Rejection')
3.2 Strawman Position Disturbance PSD
The Strawman PSD is an example of PSD enveloping. The PSD envelope is set high enough that the largest magnitude of the PSD across any frequency range is contained under this envelope. As a result of the enveloping the PSD is conservative and typically defined with only a few points.
Let the Strawman PSD structure be as follows
>> StrawmanPSD.mag = [1E6, 1E6, 1E-2]; >> StrawmanPSD.freq = [1E-3, 1, 1E4];
- the mag field is the PSD magnitude in and
- the freq field is the PSD frequencies in Hz.
3.2.1 Log Scale Interpolation of Coarse PSD definition
Often PSDs are provided with only a few points. This is done with convenience and as a measure of being conservative - i.e. the PSD contains all of the expected disturbances by enveloping the expected PSD.
Using the PSD with only a few points will almost always cause problems. As a result interpolation on a log scale is used to fill in PSD with lots of points.
If a PSD's frequency vector is defined by the vector x and the PSD's magnitude by y then interpolating a new vector of frequencies xi can be done using the interp1q command as follows
>> logX = log10(x); >> logY = log10(y); >> logXi = log10(xi); >> yi = 10.^interp1q(logX, logY, logXi);
- yi is the interpolated PSD magnitudes corresponding to the frequencies of xi.
3.3 Frequency Domain Calculation
Define a well spaced frequency vector
>> freqVec = logspace(-1, 4, 5000);
Calculate the disturbance rejection magnitude and PSD magnitude at each frequency
>> [DRmag] = bode(DR, freqVec*(2*pi)); >> logX = log10(StrawmanPSD.freq); >> logY = log10(StrawmanPSD.mag); >> logXi = log10(freqVec); >> yi = 10.^interp1q(logX(:), logY(:), logXi(:));
Note that the bode command requires the units be but the PSD is defined in Hz. As a result freqVec is converted (using 2*pi) for the bode command.
Calculate the steady-state response and plot
>> resp = yi .* DRmag(:).^2; >> figure >> semilogx(freqVec, yi, 'b', freqVec, resp, 'r-.') >> loglog(freqVec, yi, 'b', freqVec, resp, 'r-.') >> legend('Input', 'Response', 'Location', 'NorthEast') >> grid on >> ylabel('Magnitude (units^2/Hz)') >> xlabel('Frequency (Hz)') >> title('Frequency Domain Results for Strawman PSD')
Figure 5 has the resulting output PSD compared to the input PSD.
Note that this input PSD and disturbance rejection are for illustration only. The closed loop transfer function was formed from a simple 2nd order equation with little peaking. The rising slope of the disturbance rejection matches the falling slope of the input PSD leading to a flat spot in the steady state response. Real systems will have perculiarites and slopes probably won't match up so the resulting PSD for the steady state response will look different from this one.
- Sun, Jian-Qiao (2006). Stochastic Dynamics and Control, Volume 4. Amsterdam: Elsevier Science. ISBN 0444522301.