Repository: python-control/python-control

Issue:
stability_margins doesn't find gain=1 crossover frequency correctly for DT systems
stability_margins doesn't find gain=1 crossover frequency correctly for DT systems, but it seems to mostly do so for CT systems. I think the probelm may be in how `margins.poly_z_mag1_crossing` is finding roots, but there were a few  steps in that function I couldn't follow, e.g. why H(z)*H(1/z)=1 implies |H(z)| = 1. 

Example: 
```python
import control as ct
omegan = 2*np.pi
zeta = 0.2
plant = 1.1*ct.tf(omegan**2, [1, 2*zeta*omegan, omegan**2])
Ts = .05
plantdisc = ct.c2d(plant,Ts,'zoh')
#_ = ct.bode_plot(plant, margins=True) # gives correct gain=1 crossover freq
_ = ct.bode_plot(plantdisc, margins=True) # no gain=1/no phase margin given
```
<img width="407" alt="image" src="https://user-images.githubusercontent.com/58706249/106208565-6f86b100-6178-11eb-874a-809bd05a899b.png">

As an example of what it should look like, here is Matlab:
```matlab
omegan = 2*pi;
zeta = 0.2;
plant = 1.1*tf(omegan^2, [1, 2*zeta*omegan, omegan^2])
Ts = .05;
plantdisc = c2d(plant,Ts,'zoh')
margin(plantdisc)
```
<img width="476" alt="image" src="https://user-images.githubusercontent.com/58706249/106208756-bb395a80-6178-11eb-8bb2-95ee6d928a3d.png">

Discussion:

Patch:
diff --git a/control/margins.py b/control/margins.py
index 20da2a87..fc2a20d7 100644
--- a/control/margins.py
+++ b/control/margins.py
@@ -156,13 +156,14 @@ def _poly_z_real_crossing(num, den, num_inv_zp, den_inv_zq, p_q, dt, epsw):
     return z, w
 
 
-def _poly_z_mag1_crossing(num, den, num_inv, den_inv, p_q, dt, epsw):
+def _poly_z_mag1_crossing(num, den, num_inv_zp, den_inv_zq, p_q, dt, epsw):
     # |H(z)| = 1, H(z)*H(1/z)=1, num(z)*num(1/z) == den(z)*den(1/z)
-    p1 = np.polymul(num, num_inv)
-    p2 = np.polymul(den, den_inv)
+    p1 = np.polymul(num, num_inv_zp)
+    p2 = np.polymul(den, den_inv_zq)
     if p_q < 0:
+        # * z**(-p_q)
         x = [1] + [0] * (-p_q)
-        p2 = np.polymul(p2, x)
+        p1 = np.polymul(p1, x)
     z = np.roots(np.polysub(p1, p2))
     eps = np.finfo(float).eps**(1 / len(p2))
     z, w = _z_filter(z, dt, eps)
@@ -171,7 +172,7 @@ def _poly_z_mag1_crossing(num, den, num_inv, den_inv, p_q, dt, epsw):
     return z, w
 
 
-def _poly_z_wstab(num, den, num_inv, den_inv, p_q, dt, epsw):
+def _poly_z_wstab(num, den, num_inv_zp, den_inv_zq, p_q, dt, epsw):
     # Stability margin: Minimum distance to -1
 
     # TODO: Find a way to solve for z or omega analytically with given


Test patch:
diff --git a/control/tests/margin_test.py b/control/tests/margin_test.py
index 4965bbe8..fbd79c60 100644
--- a/control/tests/margin_test.py
+++ b/control/tests/margin_test.py
@@ -338,13 +338,19 @@ def test_zmore_stability_margins(tsys_zmore):
     'cnum, cden, dt,'
     'ref,'
     'rtol',
-    [([2], [1, 3, 2, 0], 1e-2,  # gh-465
-      (2.9558, 32.8170, 0.43584, 1.4037, 0.74953, 0.97079),
-      0.1  # very crude tolerance, because the gradients are not great
-      ),
-     ([2], [1, 3, 3, 1], .1,  # 2/(s+1)**3
-      [3.4927, 69.9996, 0.5763, 1.6283, 0.7631, 1.2019],
-      1e-3)])
+    [( # gh-465
+      [2], [1, 3, 2, 0], 1e-2,  
+      [2.9558, 32.390, 0.43584, 1.4037, 0.74951, 0.97079],
+      2e-3), # the gradient of the function reduces numerical precision
+     ( # 2/(s+1)**3
+      [2], [1, 3, 3, 1], .1,  
+      [3.4927, 65.4212, 0.5763, 1.6283, 0.76625, 1.2019],
+      1e-4),
+     ( # gh-523
+      [1.1 * 4 * np.pi**2], [1, 2 * 0.2 * 2 * np.pi,  4 * np.pi**2], .05,
+      [2.3842, 18.161, 0.26953, 11.712, 8.7478, 9.1504],
+      1e-4),
+     ])
 def test_stability_margins_discrete(cnum, cden, dt, ref, rtol):
     """Test stability_margins with discrete TF input"""
     tf = TransferFunction(cnum, cden).sample(dt)