# How to simulate control loop with sustained oscillation

8 views (last 30 days)
banghyun jo on 21 Mar 2024
Edited: Paul on 18 Apr 2024
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.

Balavignesh on 17 Apr 2024
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.
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 on 18 Apr 2024
Edited: Paul on 18 Apr 2024
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');
set_param(hconstant,'Value','0');
set_param(hgain,'Gain','b');
set_param(hplant,'Zeros','Gs.Z{1}');
set_param(hplant,'Poles','Gs.P{1}');
set_param(hplant,'Gain','Gs.K');
set_param(hout,'VariableName','y');
PH = 'PortHandles';
get_param(hconstant,PH).Outport,...
get_param(herror,PH).Inport(1));
get_param(hplant,PH).Outport,...
get_param(herror,PH).Inport(2));
get_param(herror,PH).Outport,...
get_param(hsaturation,PH).Inport(1));
get_param(hsaturation,PH).Outport,...
get_param(hgain,PH).Inport(1));
get_param(hgain,PH).Outport,...
get_param(hic,PH).Inport(1));
get_param(hic,PH).Outport,...
get_param(hplant,PH).Inport(1));
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.