3
$\begingroup$

I have the following Linear State Space Model. It is pasted in Normal form for it to be readable.

  provalin={{{0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {-355743., -2708.82, 
   4.38556*10^6, -23134.7, 3.24779*10^7, 252220., 1.9571*10^8, 
   1.54331*10^6, 1.3344*10^9, 1.04062*10^7, 195127., -60114.7}, {0, 0,
    0, 1, 0, 0, 0, 0, 0, 0, 0, 
   0}, {-1740.54, -2.07579, -37052.7, -188.127, 124023., 961.386, 
   2.47582*10^6, 19977.9, 1.53934*10^7, 
   121345., -5620.58, -784.516}, {0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 
   0}, {206.237, 0.144667, 0, 0, -20690.5, -378.64, 0, 0, 0, 0, 0, 
   0}, {0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0}, {0, 0, 291.405, 0.45, 0, 
   0, -18814.8, -359.96, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 
   0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, -20110., -372.677, 0, -1.}, {0, 0, 
   0, 0, 0, 0, 0, 0, 0, 0, 0, 1}, {100.113, 0.0743233, 3035.71, 
   7.23123, 16113., 124.954, 272165., 2192.23, 1.70504*10^6, 
   13428.4, -547.949, -86.0336}}, {{0}, {449.745}, {0}, {-0.412829}, \
{0}, {0}, {0}, {0}, {0}, {0}, {0}, {0.00736309}}, {{0, 0, 0, 0, 0, 0, 
   0, 0, 0, 0, 1, 0}}, {{0}}};

I've tried to construct the PSD with the following code.

   slv = OutputResponse[
   StateSpaceModel[provalin], {0.1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {t, 0, 200}];
listat = Flatten[Table[slv /. t -> tt, {tt, 0, 200, 0.1}]];
Periodogram[listat, ScalingFunctions -> "Absolute"]

I've obtained this way the following plot:

enter image description here

Though I expected to have only peaks without that initial slope. How is it? Am I wrongly setting the PSD? What are the right strategies to get this done cleanly?

$\endgroup$
4
  • $\begingroup$ Could you post the code you used to make the StateSpaceModel? $\endgroup$ – Nasser Jan 12 '17 at 23:05
  • $\begingroup$ It is terribly long and uselessly complicated. The above Linear model is a linearization of a 12 page long nonlinear system with quite complicated parameter dependencies. Have you any problems with the Normal form @Nasser? $\endgroup$ – Mirko Aveta Jan 13 '17 at 7:43
  • $\begingroup$ I just wanted to see if you used sampling rate there or not, that is all. But I have not used Periodogram before, so not sure now what the issue is. Hopefully someone else will know. $\endgroup$ – Nasser Jan 13 '17 at 7:52
  • $\begingroup$ @Nasser I've sampled directly the output response of the system (see 'listat' in my code). If you have any alternatives to Periodogram I'll be happy anyway! :) $\endgroup$ – Mirko Aveta Jan 13 '17 at 7:53
2
+50
$\begingroup$

I think the problem is that you have not computed the PSD from a steady state.

sys = StateSpaceModel[provalin];
naturalFreq = Abs[TransferFunctionPoles[sys]]/(2 Pi) // Flatten
slv = Chop[OutputResponse[sys, 0.1*UnitStep[t], t], 10^-19];
p1 = Plot[slv, {t, 0, 20}, PlotRange -> Automatic, 
          ImageSize -> Medium, PlotLabel -> "not steady state"];
p2 = Plot[slv, {t, 10, 20}, PlotRange -> Automatic, 
          ImageSize -> Medium, PlotLabel -> "steady state"];
GraphicsGrid[{{p1, p2}}]

enter image description here

Td = 0.001;
listat = Flatten[Table[slv /. t -> tt, {tt, 100, 200, Td}]];
Periodogram[listat, PlotRange -> All, SampleRate -> 1/Td, 
Epilog -> {Green, Dashed, 
    Line[{{naturalFreq[[-1]], -150}, {naturalFreq[[-1]], 0}}]}, 
    ImageSize -> Large]

enter image description here

$\endgroup$
2
  • $\begingroup$ May I ask why you divided the natural frequency by 2Pi? $\endgroup$ – Mirko Aveta Jan 23 '17 at 18:52
  • 1
    $\begingroup$ The absolute value of the pole s is the natural frequency in rad/s $\endgroup$ – OhmSweetOhm Jan 26 '17 at 20:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.