How to simulate control loop with sustained oscillation

Question

banghyun jo 2024-3-21

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2097096-how-to-simulate-control-loop-with-sustained-oscillation

编辑： Paul 2024-4-18

在 MATLAB Online 中打开

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.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Balavignesh 2024-4-17

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2097096-how-to-simulate-control-loop-with-sustained-oscillation#answer_1442751

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:

'nyquistplot' function : https://www.mathworks.com/help/ident/ref/dynamicsystem.nyquistplot.html
'freqresp' function : https://www.mathworks.com/help/rf/ref/freqresp.html

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Answer 2

Paul 2024-4-18

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2097096-how-to-simulate-control-loop-with-sustained-oscillation#answer_1443471

编辑：Paul 2024-4-18

在 MATLAB Online 中打开

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');