ninteger.m

国外人编写的一个分数阶系统的控制工具箱

            ylabel('Step response')
        case 6 % pole and zero placement
            [zeros_open_loop, poles_open_loop] = zpkdata(handles.open_loop, 'v');
            axes(handles.plot)
            plot(real(zeros_open_loop), imag(zeros_open_loop), 'ob',...
                real(poles_open_loop), imag(poles_open_loop), 'xr')
            if handles.open_loop.Ts == 0
                grid
            else
                zgrid
            end
            xlabel('Real')
            ylabel('Imaginary')
    end
else % plant parameter variations were given
    Ts = handles.open_loop.Ts;
    variable = handles.open_loop.variable;
    switch plot_type
        case 1 % Bode diagram
            w = logspace(log10(wl), log10(wh));
            [gain, phase] = bode(handles.open_loop, w);
            gain = 20*log10(squeeze(gain));
            phase = squeeze(phase);
            % gain variations
            [gain2, phase2] =...
                bode(handles.open_loop * (k+vargain(1)) / k, w);
            gain2 = 20*log10(squeeze(gain2));
            phase2 = squeeze(phase2);
            [temp1, temp2] =...
                bode(handles.open_loop * (k+vargain(2)) / k, w);
            gain2 = [gain2, 20*log10(squeeze(temp1))];
            phase2 = [phase2, squeeze(temp2)];
            % zero variations
            for i=1:length(zeros)
                if imag(zeros(i)) > 0 % this is a complex zero that was already cared for
                    continue % skip this zero
                end
                for j=1:2
                    if imag(zeros(i))==0
                        correction = zpk(zeros(i)+varzeros(i,j), zeros(i), 1,...
                            Ts, 'Variable', variable);
                    else
                        correction = zpk([zeros(i)+varzeros(i,j)...
                                conj(zeros(i))+conj(varzeros(i,j))],...
                            [zeros(i) conj(zeros(i))], 1, Ts, 'Variable', variable);
                    end
                    [temp1, temp2] =...
                        bode(tf(ss(handles.open_loop * correction)), w);
                    gain2 = [gain2, 20*log10(squeeze(temp1))];
                    phase2 = [phase2, squeeze(temp2)];
                end
            end
            % pole variations
            for i=1:length(poles)
                if imag(poles(i)) > 0 % this is a complex zero that was already cared for
                    continue % skip this zero
                end
                for j=1:2
                    if imag(poles(i))==0
                        correction = zpk(poles(i)+varpoles(i,j), poles(i), 1,...
                            Ts, 'Variable', variable);
                    else
                        correction = zpk([poles(i)+varpoles(i,j)...
                                conj(poles(i))+conj(varpoles(i,j))],...
                            [poles(i) conj(poles(i))], 1, Ts, 'Variable', variable);
                    end
                    [temp1, temp2] =...
                        bode(tf(ss(handles.open_loop * correction)), w);
                    gain2 = [gain2, 20*log10(squeeze(temp1))];
                    phase2 = [phase2, squeeze(temp2)];
                end
            end
            % plot
            axes(handles.plot1)
            semilogx(w, gain2, 'm-', w, gain, 'b-')
            grid
            xlabel('Frequency (rad/s)')
            ylabel('Gain (dB)')
            axes(handles.plot2)
            semilogx(w, phase2, 'm-', w, phase, 'b-')
            grid
            ylabel('Phase (degrees)')
        case 2 % Nichols diagram
            w = logspace(log10(wl), log10(wh));
            [gain, phase] = bode(handles.open_loop, w);
            gain = 20*log10(squeeze(gain));
            phase = squeeze(phase);
            % gain variations
            [gain2, phase2] =...
                bode(handles.open_loop * (k+vargain(1)) / k, w);
            gain2 = 20*log10(squeeze(gain2));
            phase2 = squeeze(phase2);
            [temp1, temp2] =...
                bode(handles.open_loop * (k+vargain(2)) / k, w);
            gain2 = [gain2, 20*log10(squeeze(temp1))];
            phase2 = [phase2, squeeze(temp2)];
            % zero variations
            for i=1:length(zeros)
                if imag(zeros(i)) > 0 % this is a complex zero that was already cared for
                    continue % skip this zero
                end
                for j=1:2
                    if imag(zeros(i))==0
                        correction = zpk(zeros(i)+varzeros(i,j), zeros(i), 1,...
                            Ts, 'Variable', variable);
                    else
                        correction = zpk([zeros(i)+varzeros(i,j)...
                                conj(zeros(i))+conj(varzeros(i,j))],...
                            [zeros(i) conj(zeros(i))], 1, Ts, 'Variable', variable);
                    end
                    [temp1, temp2] =...
                        bode(tf(ss(handles.open_loop * correction)), w);
                    gain2 = [gain2, 20*log10(squeeze(temp1))];
                    phase2 = [phase2, squeeze(temp2)];
                end
            end
            % pole variations
            for i=1:length(poles)
                if imag(poles(i)) > 0 % this is a complex zero that was already cared for
                    continue % skip this zero
                end
                for j=1:2
                    if imag(poles(i))==0
                        correction = zpk(poles(i)+varpoles(i,j), poles(i), 1,...
                            Ts, 'Variable', variable);
                    else
                        correction = zpk([poles(i)+varpoles(i,j)...
                                conj(poles(i))+conj(varpoles(i,j))],...
                            [poles(i) conj(poles(i))], 1, Ts, 'Variable', variable);
                    end
                    [temp1, temp2] =...
                        bode(tf(ss(handles.open_loop * correction)), w);
                    gain2 = [gain2, 20*log10(squeeze(temp1))];
                    phase2 = [phase2, squeeze(temp2)];
                end
            end
            % plot
            axes(handles.plot)
            plot(phase2, gain2, 'm-', phase, gain, 'b-')
            ngrid
            xlabel('Phase (degrees)')
            ylabel('Gain (dB)')
        case 3 % Nyquist diagram
            w = logspace(log10(wl), log10(wh));
            response = squeeze(freqresp(handles.open_loop, w));
            % gain variations
            response2 =...
                squeeze(freqresp(handles.open_loop * (k+vargain(1)) / k, w));
            response2 = [response2, squeeze(freqresp(handles.open_loop *...
                        (k+vargain(2)) / k, w))];
            % zero variations
            for i=1:length(zeros)
                if imag(zeros(i)) > 0 % this is a complex zero that was already cared for
                    continue % skip this zero
                end
                for j=1:2
                    if imag(zeros(i))==0
                        correction = zpk(zeros(i)+varzeros(i,j), zeros(i), 1,...
                            Ts, 'Variable', variable);
                    else
                        correction = zpk([zeros(i)+varzeros(i,j)...
                                conj(zeros(i))+conj(varzeros(i,j))],...
                            [zeros(i) conj(zeros(i))], 1, Ts, 'Variable', variable);
                    end
                    response2 = [response2,...
                            squeeze(freqresp(handles.open_loop * correction,...
                            w))];
                end
            end
            % pole variations
            for i=1:length(poles)
                if imag(poles(i)) > 0 % this is a complex zero that was already cared for
                    continue % skip this zero
                end
                for j=1:2
                    if imag(poles(i))==0
                        correction = zpk(poles(i)+varpoles(i,j), poles(i), 1,...
                            Ts, 'Variable', variable);
                    else
                        correction = zpk([poles(i)+varpoles(i,j)...
                                conj(poles(i))+conj(varpoles(i,j))],...
                            [poles(i) conj(poles(i))], 1, Ts, 'Variable', variable);
                    end
                    response2 = [response2,...
                            squeeze(freqresp(handles.open_loop * correction,...
                            w))];
                end
            end
            % plot
            axes(handles.plot)
            plot(real(response2), imag(response2), 'm-',...
                real(response), imag(response), 'b-', -1, 0, 'r.')
            grid
            xlabel('Real')
            ylabel('Imaginary')
        case 4 % impulse response
            if Ts>0
                t = 0:Ts:tfinal;
            else
                t = 0:tfinal/100:tfinal;
            end
            out = impulse(handles.open_loop, t);
            % gain variations
            out2 = impulse(handles.open_loop * (k+vargain(1)) / k, t);
            out2 = [out2,...
                    impulse(handles.open_loop * (k+vargain(2)) / k, t)];
            % zero variations
            for i=1:length(zeros)
                if imag(zeros(i)) > 0 % this is a complex zero that was already cared for
                    continue % skip this zero
                end
                for j=1:2
                    if imag(zeros(i))==0
                        correction = zpk(zeros(i)+varzeros(i,j), zeros(i), 1,...
                            Ts, 'Variable', variable);
                    else
                        correction = zpk([zeros(i)+varzeros(i,j)...
                                conj(zeros(i))+conj(varzeros(i,j))],...
                            [zeros(i) conj(zeros(i))], 1, Ts, 'Variable', variable);
                    end
                    out2 = [out2,...
                            impulse(handles.open_loop * correction, t)];
                end
            end
            % pole variations
            for i=1:length(poles)
                if imag(poles(i)) > 0 % this is a complex zero that was already cared for
                    continue % skip this zero
                end
                for j=1:2
                    if imag(poles(i))==0
                        correction = zpk(poles(i)+varpoles(i,j), poles(i), 1,...
                            Ts, 'Variable', variable);
                    else
                        correction = zpk([poles(i)+varpoles(i,j)...
                                conj(poles(i))+conj(varpoles(i,j))],...
                            [poles(i) conj(poles(i))], 1, Ts, 'Variable', variable);
                    end
                    out2 = [out2,...
                            impulse(handles.open_loop * correction, t)];
                end
            end
            % plot
            axes(handles.plot)
            plot(t, out2, 'm-', t, out, 'b-')
            grid
            xlabel('Time (s)')
            ylabel('Impulse response')
        case 5 % step response
            if Ts>0
                t = 0:Ts:tfinal;
            else
                t = 0:tfinal/100:tfinal;
            end
            out = step(handles.open_loop, t);
            % gain variations
            out2 = step(handles.open_loop * (k+vargain(1)) / k, t);
            out2 = [out2,...
                    step(handles.open_loop * (k+vargain(2)) / k, t)];
            % zero variations
            for i=1:length(zeros)
                if imag(zeros(i)) > 0 % this is a complex zero that was already cared for
                    continue % skip this zero
                end
                for j=1:2
                    if imag(zeros(i))==0
                        correction = zpk(zeros(i)+varzeros(i,j), zeros(i), 1,...
                            Ts, 'Variable', variable);
                    else
                        correction = zpk([zeros(i)+varzeros(i,j)...
                                conj(zeros(i))+conj(varzeros(i,j))],...
                            [zeros(i) conj(zeros(i))], 1, Ts, 'Variable', variable);
                    end
                    out2 = [out2,...
                            step(handles.open_loop * correction, t)];
                end
            end
            % pole variations
            for i=1:length(poles)
                if imag(poles(i)) > 0 % this is a complex zero that was already cared for
                    continue % skip this zero
                end
                for j=1:2
                    if imag(poles(i))==0
                        correction = zpk(poles(i)+varpoles(i,j), poles(i), 1,...
                            Ts, 'Variable', variable);
                    else
                        correction = zpk([poles(i)+varpoles(i,j)...
                                conj(poles(i))+conj(varpoles(i,j))],...
                            [poles(i) conj(poles(i))], 1, Ts, 'Variable', variable);
                    end
                    out2 = [out2,...
                            step(handles.open_loop * correction, t)];
                end
            end
            % plot
            axes(handles.plot)
            plot(t, out2, 'm-', t, out, 'b-')
            grid
            xlabel('Time (s)')
            ylabel('Step response')
        case 6 % pole and zero placement
            [zeros_open_loop, poles_open_loop] = zpkdata(handles.open_loop, 'v');
            % zero variations
            zeros2 = [];
            for i=1:length(zeros)
                if imag(zeros(i))>0
                    zeros2 = [zeros2; zeros(i)+varzeros(i-1,1);...
                            zeros(i)+varzeros(i-1,2)];
                else
                    zeros2 = [zeros2; zeros(i)+varzeros(i,1);...
                            zeros(i)+varzeros(i,2)];
                end
            end
            % pole variations
            poles2 = [];
            for i=1:length(poles)
                if imag(poles(i))>0
                    poles2 = [poles2; poles(i)+conj(varpoles(i-1,1));...
                            poles(i)+conj(varpoles(i-1,2))];
                else
                    poles2 = [poles2; poles(i)+varpoles(i,1);...
                            poles(i)+varpoles(i,2)];
                end
            end
            % plot
            axes(handles.plot)
            plot(real(zeros2), imag(zeros2), 'om',...
                real(poles2), imag(poles2), 'xg',...
                real(zeros_open_loop), imag(zeros_open_loop), 'ob',...
                real(poles_open_loop), imag(poles_open_loop), 'xr')
            if handles.open_loop.Ts == 0
                grid
            else
                zgrid
            end
            xlabel('Real')
            ylabel('Imaginary')
    end
end