How to simulate control loop with sustained oscillation

21 次查看（过去 30 天）
banghyun jo
banghyun jo 2024-3-21
编辑： Paul 2024-4-18
Ran in:
Hello.
I want to simulate below example for limit cycle
Gs = zpk([], [0 -1 -9], [9]) % open-loop system
Gs = 9 ------------- s (s+1) (s+9) Continuous-time zero/pole/gain model.
sys = feedback(Gs, 1) % closed-loop system
sys = 9 ---------------------------------- (s+9.121) (s^2 + 0.8785s + 0.9867) Continuous-time zero/pole/gain model.
figure(1)
nyquist(sys); axis([-0.4 0.1 -0.05 0.04]); grid on
how can i plot only positive value of frequency and nonlinear locus(-1/N(A))?
I find frequency(w) =3, amplitud(A) = 2/(5*pi)
so, I want to simulate u(t) = A*sin(w*t)
len = 101; % desired 101 elements in the time vector
Time= 10;
t = linspace(0,Time,len); % time vector from 0 to 10 seconds
w = 3; % frequency of the sine wave (in Hertz, not rad/s)
A = 2/(5*pi); % amplitude of the sine wave
u = A*sin(w*t); % sinusoidal input to closed-loop system
figure(2)
lsim(sys, u, t), grid on
was it right result? I'm not sure about the simulation with u(t) and feedback loop.

请先登录，再回答此问题。

回答（2 个）

Balavignesh
Balavignesh 2024-4-17
Hi Banghyun,
Your approach in simulating the system with a limit cycle using MATLAB looks almost correct to me.
The commant 'nyquist(sys); axis([-0.4 0.1 -0.05 0.04]); grid on' correctly plots the Nyquist diagram of the closed-loop system. However, Nyquist plots inherently include both positive and negative frequencies since they show the frequency response of the system. 'nyquist' function doesn't directly support plotting only the positive frequency range or a non-linear focus like '(-1/N(A))' within its standard options. To focus on positive frequencies and specific loci, you might need to manually process and plot the frequency response data using 'nyquistplot' or 'freqresp' functions and then manipulate the plot data.
Please find the documentation links below:
To ensure the simulation results align with expectations, check the time response plot generated by 'lsim' for expected behaviors, such as periodic steady-state oscillations that would indicate a limit cycle and compare the amplitude and frequency of the oscillations in the output with the input signal to see if they match the expected limit cycle characteristics.
Hope that helps!
Balavignesh
Paul
Paul 2024-4-18
编辑：Paul 2024-4-18
Ran in:
Hi banghyun jo,
Describing function analysis relies on the Nyquist plot of the open-loop system, not the closed-loop system. If -1/N(A) crosses the Nyquist plot of the open-loop system, then we have the potential for sustained limit cycle.
Define the plant
Gs = zpk([], [0 -1 -9], [9]);
Define and plot the describing function of the ideal relay nonlinear element wtth b = 1
N = @(A,b) 4*b./pi./A;
A = 0:.01:.5;b=1;
figure
plot(A,N(A,b)),grid
xlabel('A')
ylabel('N(A)')
I'll add a line at N(A) = 10. We see that N(A) = 10 corresponds to A = 0.125 (close enough), which will be important later.
yline(10)
Plot the Nyquist plot of the open-loop system and overlay with -1/N(A).
figure
nyquist(Gs)
hold on
plot(real(-1./N(A,b)),0*A,'g')
ylim([-1 1])
-1/N(A) traverses from the origin to the left and crosses the Nyquist plot of G(s). A bit of analysis would show that we should expect a sustained limit cycle.
-1/N(A) cross the Nyquist plot on the real axis where the phase of Gs(s) is -180 deg. Plot the Bode plot of Gs to determine the frequency and gain at that point.
figure
margin(Gs),grid
The frequency is w = 3 rad/sec and the gain is -20 dB where the phase crosses -180 deg. For a neutrally stable system, i.e., sustained limit cycle, the effective gain of the nonlinear element must be 20 dB, or 10 in absolute terms. Recall from above that N(0.125) = 10, hence the expected amplitude of the oscillation is
Aexpected = 0.125;
The frequency (Hz) of the oscillation is determined from the 180 deg crossover frequency
Fexpected = 3/2/pi;
Therefore the period of the oscillation should be
Texpected = 1./Fexpected
Texpected = 2.0944
We might be able to simulate this system with one of the ode solvers, but perhaps it's better to use Simulink and take advantage of its zero crossing detection cabability. The following code builds a model of your block diagram (and includes an initial condition on the plant input to stimulate the system).
bdclose('all');
hsys = new_system('dfnctest');
hconstant = add_block('simulink/Sources/Constant','dfnctest/ref');
set_param(hconstant,'Value','0');
herror = add_block('simulink/Math Operations/Subtract','dfnctest/error');
hsaturation = add_block('simulink/Math Operations/Sign','dfnctest/sign');
hgain = add_block('simulink/Math Operations/Gain','dfnctest/gain');
set_param(hgain,'Gain','b');
hic = add_block('simulink/Signal Attributes/IC','dfnctest/IC');
hplant = add_block('simulink/Continuous/Zero-Pole','dfnctest/plant');
set_param(hplant,'Zeros','Gs.Z{1}');
set_param(hplant,'Poles','Gs.P{1}');
set_param(hplant,'Gain','Gs.K');
hout = add_block('simulink/Sinks/To Workspace','dfnctest/output');
set_param(hout,'VariableName','y');
PH = 'PortHandles';
add_line(hsys, ...
get_param(hconstant,PH).Outport,...
get_param(herror,PH).Inport(1));
add_line(hsys, ...
get_param(hplant,PH).Outport,...
get_param(herror,PH).Inport(2));
add_line(hsys, ...
get_param(herror,PH).Outport,...
get_param(hsaturation,PH).Inport(1));
add_line(hsys, ...
get_param(hsaturation,PH).Outport,...
get_param(hgain,PH).Inport(1));
add_line(hsys, ...
get_param(hgain,PH).Outport,...
get_param(hic,PH).Inport(1));
add_line(hsys, ...
get_param(hic,PH).Outport,...
get_param(hplant,PH).Inport(1));
add_line(hsys, ...
get_param(hplant,PH).Outport,...
get_param(hout,PH).Inport(1));
% open_system(hsys) % to open on the Simulink canvas
cset = getActiveConfigSet(hsys);
set_param(cset,'StopTime','20');
set_param(cset,'Solver','VariableStepAuto');
set_param(cset,'MaxStep','0.01');
Simulate the model and plot the output
out = sim('dfnctest',cset);
figure
plot(out.y),grid
As predicted, a sustained oscillation with period of 2.1 sec and amplitude of 0.125.

类别

在 Help Center 和 File Exchange 中查找有关 Time and Frequency Domain Analysis 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by