1 function nyqlog(sys) 2 3 % Nyqlog makes a polar plot of the open 4 % loop transfer function h0(s) with 5 % |h0(s)| on a dB scale. It supports 6 % only continuous and monovariable systems. 7 8 % Original MATLAB version: 9 % 10 % Trondheim, June 2009 11 % Trond Andresen <Trond.Andresen@itk.ntnu.no> 12 % Department of Engineering Cybernetics, 13 % The Norwegian University of Science and Technology, 14 % N-7034 Trondheim, NORWAY. 15 16 % Port to Octave: 17 % 18 % Buenos Aires, May 2026 19 % Martin Klöckner <mklockner@fi.uba.ar> 20 % Facultad de Ingenieria, 21 % Universidad de Buenos Aires, 22 % Argentina. 23 24 % May be distributed freely for non-commercial use, 25 % but please leave the above info unchanged, for 26 % credit and feedback purposes. 27 28 % *********************************************** 29 % System examples for copying and pasting on the 30 % command line: 31 % ------------------- 32 % sys=tf(1, conv([1 0 0], [1 1])); 33 % sys=tf(conv([5 1],[1 1]), [1 0 0 0],'ioDelay',1); 34 % sys=tf(1,conv([1 -0.01 0 1 3 -0.1 7],[1 -0.05 0.6])); 35 % sys=tf(1,[1 0.00000001 1+0.00000001^2]); 36 % sys=tf(1,[1 0]); 37 % w0=10; zeta=0.1; sys=tf([1 1],[1/w0^2 2*zeta/w0 1]); 38 % w0=5; zeta=0.1; sys=tf(1,[1/w0^2 2*zeta/w0 1]); 39 % sys=0.1*tf([10 1], [5 6 1 0 0]); 40 % sys=1*tf([1],conv([1 -0.1],[1 0])); 41 % sys=tf([1 0 25],conv([1 0 1 0 0],[1 0 4])); 42 % sys=tf([1 0 25],[1 0 1 ]); 43 % sys=tf([1 0 1], [1 0 0 0]); 44 % sys = zpk([-1/3 -1/2], [0 -0.02 -0.1 -2 -10], 48) 45 46 % Table 9.6 - 3: 47 % sys=10*tf(1,conv(conv([1000 1],[10 1]),[1 1])) 48 % Table 9.6 - 9: 49 % sys=tf(0.001,[50 1 0 0]) 50 % Table 9.6 - 12: 51 % sys=0.05*tf([5 1], [1 0 0 0]) 52 % Table 9.6 - 13: 53 % sys=tf(conv([5 1],[1 1]), [1 0 0 0]); 54 % Table 9.6 - 14: 55 % sys=200*tf(conv([3 1],[2 1]),conv(conv([50 1 0],[10 1]),conv([0.5 2],[0.1 1]))) 56 % Table 9.6 - 15: 57 % sys=10*tf([25 1],conv([1 2 0 0],[1 1])) 58 59 % Ex. 8.15 in REGULERINGSTEKNIKK (in Norwegian) by Balchen, Andresen, Foss: 60 % a=0.2; T =1; Kp= 1; sys=tf(Kp,conv([1 -a],[T 1])); 61 % ************************************************* 62 63 % Checking system order: 64 [num, den] = tfdata(sys, "vector"); 65 66 idx = find(num != 0, 1); 67 if !isempty(idx) 68 num = num(idx:end); 69 endif 70 71 if !issiso(sys) 72 error("Only monovariable systems allowed"); 73 endif 74 75 % Checking for delay: 76 try 77 dly = get(sys, "ioDelay"); 78 catch 79 dly = 0; 80 end_try_catch 81 82 % Checking that the system is not discrete: 83 Ts = get(sys, "tsam"); 84 85 if isempty(Ts) || Ts > 0 86 error("Only continuous systems allowed"); 87 endif 88 89 90 % Charting poles and zeroes, system dimension; 91 % Sorting poles by Im-value in ascending order: 92 93 sysdim = size(den); 94 sysdim = sysdim(2)-1; 95 if (sysdim == 0) 96 error('Denominator order of zero not allowed in Nyqlog'); 97 end 98 99 numdim = size(num); 100 numdim = numdim(2)-1; 101 if (numdim > sysdim) 102 error('Denominator order must be >= nominator order'); 103 end 104 105 poles = roots(den); 106 re_poles = real(poles); 107 % This is used below to check stability of the closed-loop system 108 openloop_rhp_poles = sum(re_poles > 0); 109 impoles = imag(poles); 110 [vhlp,im_ndx] = sort(impoles); 111 poles(:) = poles(im_ndx(:)); 112 abpoles = abs(poles); 113 114 repoles = real(poles); 115 impoles = imag(poles); 116 117 % No. of poles in the origin and on the imaginary axis, if any; 118 Np_origin = 0; 119 Np_imag = 0; 120 for k = 1 : sysdim 121 if (repoles(k) == 0) 122 if (impoles(k) > 0) 123 Np_imag = Np_imag+1; 124 elseif (impoles(k) == 0) 125 Np_origin = Np_origin+1; 126 else 127 ; 128 end 129 end 130 end 131 132 if Np_imag 133 impoles = impoles(end+1-Np_imag:end); 134 end 135 136 % Sorting zeroes by Im-value in ascending order: 137 abzeroes = 0; 138 if (numdim) 139 zeroes = roots(num); 140 imzeroes = imag(zeroes); 141 [vhlp, im_ndx] = sort(imzeroes); 142 zeroes(:) = zeroes(im_ndx(:)); 143 abzeroes = abs(zeroes); 144 rezeroes = real(zeroes); 145 imzeroes = imag(zeroes); 146 end 147 148 % Finding smallest distance > 0 to pole or zero 149 % to decide radius for half circles around imaginary poles; 150 eps = 1e-3; 151 R = eps; 152 153 if Np_origin 154 R = eps^(1/Np_origin); 155 end 156 157 R0 = inf; 158 if (Np_origin && Np_imag) 159 R = mindist(0, poles, R0); 160 if (numdim) 161 R = mindist(0, zeroes, R0); 162 end 163 end 164 165 R0 = R; 166 if Np_imag 167 for k = 1 : length(impoles) 168 R = mindist(impoles(k), poles, R0); 169 end 170 end 171 172 R0=R; 173 if (numdim) 174 for k = 1 : length(impoles) 175 R = mindist(impoles(k), zeroes, R0); 176 end 177 end 178 179 frac = 0.01; 180 if Np_origin 181 frac = (0.01)^(1/Np_origin); 182 end 183 184 R=frac*R; 185 186 % Calculating a tentative vector s to be used in h0(s): 187 % Ns = number of points on s. Special measures are 188 % taken if there is a time delay in sys. 189 Ns = 1000 + 150*(dly > 0); 190 wmin = R; 191 wmax = 1e6*max(cat(1, abpoles, abzeroes)); 192 193 if (wmax == 0) 194 wmax= R*1e11; 195 end 196 197 if (dly) 198 wmax = min(wmax, 4*pi/dly); 199 end 200 201 w = log(wmin):(log(wmax) - log(wmin)) / Ns:log(wmax); 202 w = exp(w'); 203 s = j.*w; 204 205 % Splicing in possible imaginary zeroes > 0 in s: 206 n = 0; 207 for k = 1 : numdim 208 if (abs(rezeroes(k)) < 1e-12 && imzeroes(k) > R) 209 n = n + 1; 210 imonlyzeroes(n) = imzeroes(k); 211 end 212 end 213 if (exist('imonlyzeroes')) 214 w1 = cat(1, w, imonlyzeroes'); 215 [w1, im_ndx] = sort(w1, 1); 216 s = j.*w1; 217 end 218 219 % Plot real and imaginary axes 220 plot([-1.5 1.5], [0 0], 221 'Color', [0.53 0.53 0.53], 222 'LineWidth', 1.0); 223 224 hold on; 225 226 plot([0 0], [-1.5 1.5], 227 'Color', [0.53 0.53 0.53], 228 'LineWidth', 1.0); 229 230 % Plotting background diagram 231 circle(0,0.5,'-'); 232 circle(0,0.6666666667,'-.'); 233 circle(0,0.833333333,'-.'); 234 circle(0,1,'-', 'LineWidth', 1.25, 'Color', [0.55 0.55 0.55]); 235 circle(0,1.16666667,'-.'); 236 circle(0,1.33333333,'-.'); 237 circle(0,1.5,'-'); 238 phase_lines(24, 0.5, 1.5, '-.'); 239 phase_lines(8, 0.5, 1.5, '-'); 240 241 242 % Plotting main graph and its mirror image: 243 s1 = s; 244 spiralfactor = 1.5; 245 spiralfactor = 1/spiralfactor; 246 247 % set(gcf, 'Color', [0.1 0.1 0.1]); 248 set(gcf, "Color", get(0, "defaultfigurecolor")); 249 s = scurve(s1, R, spiralfactor, Np_origin, Np_imag, impoles); 250 %plot(s); hold on; 251 %break; 252 %plot([0 1],[R 1]); 253 s = conj(s); 254 [zmirr, ncount] = nygraph(sys, s, 255 '--', 256 'LineWidth', 2.0, 257 'Color', get(0, "defaultAxesXColor")); 258 259 spiralfactor = 1/spiralfactor; 260 261 s = scurve(s1, R, spiralfactor, Np_origin, Np_imag, impoles); 262 [zmain, ncount] = nygraph(sys, s, '-', 263 'LineWidth', 2.0, 264 'Color', [1.0 0.0 0.0]); 265 266 fprintf(1,'Number of poles in RHP of open-loop system: %i\n', openloop_rhp_poles); 267 268 % Count net encirclements around the point -1. 269 % An encirclement is counted as positive if the direction 270 % is clockwise. 271 [ncirc,npoles_on_im_axis] = countencirc(zmirr,zmain); 272 if npoles_on_im_axis > 0 273 fprintf(1,'%d closed-loop pole(s) on the im-axis.\n',npoles_on_im_axis); 274 fprintf(1,'Graph goes through the -1 point, so\n'); 275 fprintf(1,'encirclement counting cannot be done.\n'); 276 else 277 fprintf(1,'Number of net encirclements around the -1 point: %i\n',ncirc); 278 closedloop_rhp_poles=ncirc+openloop_rhp_poles; 279 fprintf(1,'=> Number of poles in RHP of closed-loop system: %i\n',... 280 closedloop_rhp_poles); 281 if closedloop_rhp_poles > 0 282 fprintf(1,'=> Closed-loop-system is unstable\n'); 283 else 284 fprintf(1,'and no closed-loop poles on Im-axis\n'); 285 fprintf(1,'=> Closed-loop-system is asymptotically stable\n'); 286 end 287 end 288 289 bg_color = get(gcf, "color"); % RGB triplet, values 0..1 290 inv_bg_color = 1 - bg_color; % invert color 291 292 % Point (-1, 0) 293 plot(-1, 0, 'o', 294 'markersize', 6, 295 'markeredgecolor', inv_bg_color, 296 'Color', inv_bg_color); 297 298 text(0.03,-0.05,'-120','FontSize',16); 299 text(0.4141, -0.3331,'-60','FontSize',16); 300 text(0.8, -0.64,'0 dB','FontSize',16); 301 text(1.1924, -0.9575,'+60','FontSize',16); 302 303 % Plotting directional arrows: 304 for xh = [0.65 0.45 0.05] 305 % for xh = [0.8 0.15] 306 nmid=round(xh*ncount); 307 % arrow(zmain(nmid+1), zmain(nmid),'-'); 308 arrow(zmain(nmid+1), zmain(nmid), 309 '-', 310 'LineWidth', 2.5, 311 'Color', [1.0 0.0 0.0]); 312 end 313 314 % for xh = [0.75 0.6 0.4 0.2] 315 for xh = [0.03 0.45 0.90] 316 % for xh = [0.7 0.1] 317 nmid = round(xh*ncount); 318 % arrow(zmirr(nmid), zmirr(nmid+1),'-'); 319 arrow(zmirr(nmid), zmirr(nmid + 1), 320 '-', 321 'LineWidth', 2.5, 322 'Color', get(0, "defaultAxesXColor")); 323 end 324 325 % % Contours for |N|= const. may be plotted: 326 % % for instance [6 3 1 0.5 0.25 0 -0.5 -1 -3 -6] 327 328 % nlgrid([6 3]); 329 330 scalexy = axis; scalexy(3:4) = 1.00*scalexy(3:4); 331 axis(scalexy); 332 axis equal; axis off; axis tight; hold off; 333 334 % h = title("Nyquist plot with logaritmic amplitude"); 335 336 % set(h, "units", "normalized"); 337 % p = get(h, "position"); 338 % p(2) = p(2) * 1.01; % default is near 1.0 339 % set(h, "position", p); 340 341 endfunction 342 343 %****************************************** 344 %* SUB-FUNCTIONS * 345 %****************************************** 346 347 function arrow(z2, z1, varargin) 348 349 % dz=0.12*exp(j*angle(z2-z1)); 350 dz = 0.10*exp(j*angle(z2-z1)); 351 z_arrow_end1 = z2 - dz*exp(j*pi/4); 352 z_arrow_end2 = z2 - dz*exp(-j*pi/4); 353 plot([real(z2) real(z_arrow_end1)], 354 [imag(z2) imag(z_arrow_end1)], 355 'LineWidth', 1.5, 356 'Color', [1.0 0.0 0.0], 357 varargin{:}); 358 359 plot([real(z2) real(z_arrow_end2)],... 360 [imag(z2) imag(z_arrow_end2)], 361 'LineWidth', 1.5, 362 'Color', [1.0 0.0 0.0], 363 varargin{:}); 364 365 endfunction 366 367 %*********************************************** 368 function circle(zcentre, radius, varargin) 369 370 angles=0:pi/72:2*pi; 371 circ=zcentre+radius.*(cos(angles)+j.*sin(angles)); 372 373 plot(circ, 374 'LineWidth', 1.0, 375 'Color', [0.5 0.5 0.5], 376 varargin{:}); 377 378 endfunction 379 380 %*********************************************** 381 function phase_lines(n,rstart,rend,varargin) 382 hold on; 383 angles = 0:2*pi/n:2*pi; 384 lines = ones(n,2); 385 for k = 1 : n 386 zh = cos(angles(k)) + j*sin(angles(k)); 387 lines(k,1) = rstart*zh; 388 lines(k,2) = rend*zh; 389 plot(real(lines(k,:)), imag(lines(k,:)), 390 'LineWidth', 1.0, 391 'Color', [0.5 0.5 0.5], 392 varargin{:}); 393 end 394 endfunction 395 396 %*********************************************** 397 398 function [s]= scurve(s1,R,spiralfactor, Np_origin, Np_imag,impoles) 399 a = log(spiralfactor)*2/pi; 400 % R = R/spiralfactor; 401 % Calculating first arc if pole(s) in the origin: 402 % If there are one or more pure integrators s is 403 % made to do a small arc of a log spiral 404 % into the upper right quadrant from 0 to pi/2. 405 s = s1; 406 if Np_origin 407 fi = 0 : 0.02/Np_origin : 1; 408 fi = 0.5*pi*fi; 409 sarc = R*exp((a+j)*fi); 410 % merging sarc with s: 411 x1 = imag(sarc(end)); 412 k = 1; 413 414 while (imag(s1(k)) < x1) 415 k=k+1; 416 end 417 418 s1 = cat(1, sarc.', s(k:end)); 419 s = s1; 420 end 421 422 % Calculating arcs for possible pole(s) on the im-axis: 423 % For each such possible pole we generate a vector 424 % sarc describing a log spiral around the pole 425 % from -pi/2 to pi/2 into the right half plane. 426 fi = -1:0.02:1; 427 fi = 0.5*pi*fi; 428 s1 = s; 429 for m = 1:Np_imag 430 R = R/spiralfactor; 431 sarc = R*exp((a+j)*fi); 432 x1 = impoles(m)+imag(sarc(1)); 433 k = 1; 434 435 while (imag(s1(k)) < x1) 436 k = k+1; 437 end 438 439 sarc(:) = sarc(:) + j*impoles(m); 440 s1 = cat(1, s1(1:k-1), sarc.'); 441 x1 = impoles(m) + imag(sarc(end)); 442 k = 1; 443 while (imag(s(k)) < x1) 444 k = k+1; 445 end 446 s1 = cat(1, s1, s(k:end)); 447 end 448 s = s1; 449 endfunction 450 451 %*********************************************** 452 function val = evalfr(sys, x) 453 % Local replacement for evalfr(sys, x) 454 % Works for SISO transfer functions 455 [num, den] = tfdata(sys, "vector"); 456 numv = polyval(num, x); 457 denv = polyval(den, x); 458 val = numv / denv; 459 endfunction 460 461 function [zplot, ncount] = nygraph(sys, s, varargin) 462 % Preparing logarithmic polar plot data: 463 kmax = size(s); 464 kmax = kmax(1); 465 for k = 1 : kmax 466 zh(k)=evalfr(sys, s(k)); 467 end 468 469 z = zh.'; 470 absz = abs(z) + 1e-14; 471 logabs = 20.*log10(absz); 472 473 % Avoiding plot continuing for |h0| < -120dB: 474 % the vector s is then truncated. 475 for k = 1 : kmax 476 if (logabs(k) <= -120) 477 logabs(k) = -120; 478 end 479 end 480 481 ncount = length(logabs); 482 while (logabs(ncount) <= -120) 483 ncount = ncount - 1; 484 end 485 % From now on all vectors are ncount long; ncount <= size(s). 486 487 % Plotting the two conjugate halves of the polar curve; 488 logabsplot = logabs(1:ncount) ./ 120.+1; 489 zplot = z(1:ncount) .* logabsplot ./ absz(1 : ncount); 490 491 % plot(zplot, 'LineWidth', 2.2); 492 plot(zplot, varargin{:}); 493 494 endfunction 495 496 % ************************ 497 498 function [ncirc,npoles_on_im_axis] = countencirc(zmirr,zmain) 499 500 % Counts net encirclements around the point -1. 501 % An encirclement is counted as positive if the direction 502 % is clockwise. This follows Belanger (1995): 503 % "Control Engeering", Saunders College Publishing, 504 % pp 206 - 208. 505 % 506 % Bugs fixed and improvements made in Feb. 09: 507 % The function now also counts poles on the im-axis for the 508 % closed-loop system, if any. If such poles exist, this 509 % corresponds to the graph going through -1. 510 % Encirclement counting is then impossible and is 511 % disabled. 512 % Another (small) bug fixed and improvements made in June 09 513 514 eps = 1e-6; 515 zmirr(1:end) = zmirr(end:-1:1); 516 zmirr = zmirr(2:end-1); 517 zall = [zmirr;zmain;zmirr(1)]; 518 if abs(imag(zall(1))) < eps 519 zall = [zall;zall(2)]; 520 end 521 ncirc = 0; 522 npoles_on_im_axis = 0; 523 z3 = zall(end); 524 for k = 3 : length(zall) 525 z4 = z3; 526 z1 = zall(k); 527 z2 = zall(k-1); 528 z3 = zall(k-2); 529 abz1 = abs(z1+1); 530 abz2 = abs(z2+1); 531 abz3 = abs(z3+1); 532 zre1 = real(z1); 533 zre2 = real(z2); 534 535 % Checking if graph is too close to -1: 536 dl1 = fromline2minusone(z1,z2); 537 dl2 = fromline2minusone(z2,z3); 538 dl3 = fromline2minusone(z3,z4); 539 540 closest_now = abz1 > abz2 && abz3 > abz2; 541 542 if closest_now && min([dl1 dl2 dl3]) < 1e-5 && min([abz1 abz2 abz3]) < 0.001 543 npoles_on_im_axis = npoles_on_im_axis + 1; 544 end 545 % end checking if graph is too close to -1. 546 547 % Only checking for Re axis crossings to the left 548 % of minus 0.9 to avoid unnecessary work: 549 if zre1 < -0.9 550 zim1 = imag(z1); 551 zim2 = imag(z2); 552 zim3 = imag(z3); 553 554 if zim1*zim2 < 0 555 % Interpolation to find real z value at crossing: 556 delta12 = (real(z1) - real(z2)) * abs(imag(z2)) / (abs(imag(z1)) 557 + abs(imag(z2))); 558 realcross = real(z2) + delta12; 559 end 560 if zim1 > eps && zim2 < -eps 561 if realcross < -1 562 ncirc = ncirc+1;% ncirc,z1, z2, z3 563 end 564 elseif zim1 < -eps && zim2 > eps 565 if realcross < -1 566 ncirc = ncirc-1;% ncirc,z1, z2, z3 567 end 568 elseif abs(zim2) < eps && zim1 > 0 && zim3 < 0 569 if real(z2) < -1 570 ncirc = ncirc+1;% ncirc,z1, z2, z3 571 end 572 elseif abs(zim2) < eps && zim3 > 0 && zim1 < 0 573 if real(z2) < -1 574 ncirc = ncirc-1;% % ncirc,z1, z2, z3 575 end 576 else 577 end 578 end % real(z1) < -0.99 579 end 580 endfunction 581 582 % ************************ 583 584 function [dist] = mindist(point,vector,initdist) 585 % Calculates the minimum distance from a given complex number 586 % to a set of other complex numbers: 587 mdist = initdist; 588 kmax = length(vector); 589 for k = 1 : kmax 590 d0 = abs(vector(k) - point); 591 if (d0 ~= 0) 592 mdist = min(d0, mdist); 593 end 594 end 595 dist = mdist; 596 endfunction 597 598 % ************************ 599 600 function nlgrid(absNdB) 601 % absNdB = [6 3 1 0.5 0.25 0 -0.5 -1 -3 -6] 602 absNdB = absNdB'; 603 n = length(absNdB); 604 absN = 10.^(absNdB/20); 605 radii = 1./absN; 606 nangles = 200; 607 angles = 0:pi/nangles:2*pi; 608 angles = angles'; 609 for k = 1 : n 610 circ =-1.+radii(k) .* (cos(angles) + j.*sin(angles)); 611 absc = abs(circ)+1e-14; 612 logabs = 20.*log10(absc); 613 for p = 1:nangles 614 if (logabs(p) <= -120) 615 logabs(p) = -120; 616 end 617 end 618 619 logabsplot = logabs./120.+1; 620 cplot = circ.*logabsplot./absc; 621 % plot(cplot,'k-','LineWidth',0.5); 622 plot(cplot,'-','LineWidth',1.0); 623 % plot(cplot,'-','LineWidth',1.0, 'Color', [1.0 1.0 1.0]); 624 end 625 endfunction 626 627 % ******************************** 628 629 function [distline_to_minus1] = fromline2minusone(z1,z2) 630 % Calcucates the min. distance from the point -1 to the line z1,z2 631 v = z2-z1; 632 v = imag(v) - i*real(v); 633 v = v/abs(z2-z1); 634 r = z1+1; 635 d = dot([real(v) imag(v)], [real(r) imag(r)]'); 636 distline_to_minus1 = abs(d); 637 % pointonline=-1-distline_to_minus1*v; 638 endfunction 639 640 % ********************************
