How to simulate control loop with sustained oscillation

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.

Answer 1

[Balavignesh]

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

Answer 2

[Paul]

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.