Repository navigation
Passivity test is ill-conditioned #761
Description
Whack-a-mole!
Now this cvxopt test is flaky passivity_test.py::test_ispassive_edge_cases :
https://github.com/python-control/python-control/runs/7930805491?check_suite_focus=true#step:6:3450
Details
2022-08-20T10:05:12.4606678Z =================================== FAILURES ===================================
2022-08-20T10:05:12.4607612Z _________________ test_ispassive_edge_cases[test_input2-True] __________________
2022-08-20T10:05:12.4607862Z
2022-08-20T10:05:12.4608255Z test_input = (array([[-3.e+12, 0.e+00],
2022-08-20T10:05:12.4608627Z [ 0.e+00, -2.e+12]]), array([[0],
2022-08-20T10:05:12.4609085Z [1]]), array([[-1, 2]]), array([[0.]]))
2022-08-20T10:05:12.4609312Z expected = True
2022-08-20T10:05:12.4609448Z
2022-08-20T10:05:12.4609567Z @pytest.mark.parametrize(
2022-08-20T10:05:12.4609948Z "test_input,expected",
2022-08-20T10:05:12.4610188Z [((A, B, C, D*0.0), True),
2022-08-20T10:05:12.4610415Z ((A_d, B, C, D), True),
2022-08-20T10:05:12.4610640Z ((A*1e12, B, C, D*0), True),
2022-08-20T10:05:12.4610854Z ((A, B*0, C*0, D), True),
2022-08-20T10:05:12.4611073Z ((A*0, B, C, D), True),
2022-08-20T10:05:12.4611293Z ((A*0, B*0, C*0, D*0), True)])
2022-08-20T10:05:12.4611683Z def test_ispassive_edge_cases(test_input, expected):
2022-08-20T10:05:12.4611958Z A = test_input[0]
2022-08-20T10:05:12.4612178Z B = test_input[1]
2022-08-20T10:05:12.4612610Z C = test_input[2]
2022-08-20T10:05:12.4612970Z D = test_input[3]
2022-08-20T10:05:12.4613469Z sys = ss(A, B, C, D)
2022-08-20T10:05:12.4613753Z > assert(passivity.ispassive(sys) == expected)
2022-08-20T10:05:12.4613952Z
2022-08-20T10:05:12.4614077Z control/tests/passivity_test.py:115:
2022-08-20T10:05:12.4614373Z _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
2022-08-20T10:05:12.4614742Z control/passivity.py:292: in ispassive
2022-08-20T10:05:12.4615096Z return solve_passivity_LMI(sys, rho=ofp_index, nu=ifp_index) is not None
2022-08-20T10:05:12.4615456Z control/passivity.py:177: in solve_passivity_LMI
2022-08-20T10:05:12.4615750Z sol = cvx.solvers.sdp(c, Gs=Gs, hs=hs)
2022-08-20T10:05:12.4616316Z /usr/share/miniconda3/envs/test-env/lib/python3.10/site-packages/cvxopt/coneprog.py:4126: in sdp
2022-08-20T10:05:12.4616810Z sol = conelp(c, G, h, dims, A = A, b = b, primalstart = ps, dualstart = ds, kktsolver = kktsolver, options = options)
2022-08-20T10:05:12.4617643Z /usr/share/miniconda3/envs/test-env/lib/python3.10/site-packages/cvxopt/coneprog.py:1395: in conelp
2022-08-20T10:05:12.4618038Z misc.update_scaling(W, lmbda, ds, dz)
2022-08-20T10:05:12.4618322Z _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
2022-08-20T10:05:12.4618485Z
2022-08-20T10:05:12.4618810Z W = {'beta': [], 'd': <0x1 matrix, tc='d'>, 'di': <0x1 matrix, tc='d'>, 'dnl': <0x1 matrix, tc='d'>, ...}
2022-08-20T10:05:12.4619335Z lmbda = <6x1 matrix, tc='d'>, s = <13x1 matrix, tc='d'>
2022-08-20T10:05:12.4619636Z z = <13x1 matrix, tc='d'>
2022-08-20T10:05:12.4619773Z
2022-08-20T10:05:12.4619885Z def update_scaling(W, lmbda, s, z):
2022-08-20T10:05:12.4620221Z """
2022-08-20T10:05:12.4620607Z Updates the Nesterov-Todd scaling matrix W and the scaled variable
2022-08-20T10:05:12.4621042Z lmbda so that on exit
2022-08-20T10:05:12.4621247Z
2022-08-20T10:05:12.4621547Z W * zt = W^{-T} * st = lmbda.
2022-08-20T10:05:12.4621773Z
2022-08-20T10:05:12.4622145Z On entry, the nonlinear, 'l' and 'q' components of the arguments s
2022-08-20T10:05:12.4622629Z and z contain W^{-T}*st and W*zt, i.e, the new iterates in the current
2022-08-20T10:05:12.4622920Z scaling.
2022-08-20T10:05:12.4623119Z
2022-08-20T10:05:12.4623489Z The 's' components contain the factors Ls, Lz in a factorization of
2022-08-20T10:05:12.4624030Z the new iterates in the current scaling, W^{-T}*st = Ls*Ls',
2022-08-20T10:05:12.4624341Z W*zt = Lz*Lz'.
2022-08-20T10:05:12.4624645Z """
2022-08-20T10:05:12.4624832Z
2022-08-20T10:05:12.4625008Z
2022-08-20T10:05:12.4625390Z # Nonlinear and 'l' blocks
2022-08-20T10:05:12.4625606Z #
2022-08-20T10:05:12.4626008Z # d := d .* sqrt( s ./ z )
2022-08-20T10:05:12.4626266Z # lmbda := lmbda .* sqrt(s) .* sqrt(z)
2022-08-20T10:05:12.4626501Z
2022-08-20T10:05:12.4626752Z if 'dnl' in W:
2022-08-20T10:05:12.4627026Z mnl = len(W['dnl'])
2022-08-20T10:05:12.4627256Z else:
2022-08-20T10:05:12.4627460Z mnl = 0
2022-08-20T10:05:12.4627713Z ml = len(W['d'])
2022-08-20T10:05:12.4627934Z m = mnl + ml
2022-08-20T10:05:12.4628165Z s[:m] = base.sqrt( s[:m] )
2022-08-20T10:05:12.4628397Z z[:m] = base.sqrt( z[:m] )
2022-08-20T10:05:12.4629372Z
2022-08-20T10:05:12.4629630Z # d := d .* s .* z
2022-08-20T10:05:12.4629940Z if 'dnl' in W:
2022-08-20T10:05:12.4630299Z blas.tbmv(s, W['dnl'], n = mnl, k = 0, ldA = 1)
2022-08-20T10:05:12.4630693Z blas.tbsv(z, W['dnl'], n = mnl, k = 0, ldA = 1)
2022-08-20T10:05:12.4631031Z W['dnli'][:] = W['dnl'][:] ** -1
2022-08-20T10:05:12.4631420Z blas.tbmv(s, W['d'], n = ml, k = 0, ldA = 1, offsetA = mnl)
2022-08-20T10:05:12.4631840Z blas.tbsv(z, W['d'], n = ml, k = 0, ldA = 1, offsetA = mnl)
2022-08-20T10:05:12.4632171Z W['di'][:] = W['d'][:] ** -1
2022-08-20T10:05:12.4632391Z
2022-08-20T10:05:12.4632808Z # lmbda := s .* z
2022-08-20T10:05:12.4633059Z blas.copy(s, lmbda, n = m)
2022-08-20T10:05:12.4633323Z blas.tbmv(z, lmbda, n = m, k = 0, ldA = 1)
2022-08-20T10:05:12.4633622Z
2022-08-20T10:05:12.4633806Z
2022-08-20T10:05:12.4634057Z # 'q' blocks.
2022-08-20T10:05:12.4634267Z #
2022-08-20T10:05:12.4634539Z # Let st and zt be the new variables in the old scaling:
2022-08-20T10:05:12.4634792Z #
2022-08-20T10:05:12.4635010Z # st = s_k, zt = z_k
2022-08-20T10:05:12.4635228Z #
2022-08-20T10:05:12.4635560Z # and a = sqrt(st' * J * st), b = sqrt(zt' * J * zt).
2022-08-20T10:05:12.4635811Z #
2022-08-20T10:05:12.4636209Z # 1. Compute the hyperbolic Householder transformation 2*q*q' - J
2022-08-20T10:05:12.4636519Z # that maps st/a to zt/b.
2022-08-20T10:05:12.4636747Z #
2022-08-20T10:05:12.4637058Z # c = sqrt( (1 + st'*zt/(a*b)) / 2 )
2022-08-20T10:05:12.4637420Z # q = (st/a + J*zt/b) / (2*c).
2022-08-20T10:05:12.4637667Z #
2022-08-20T10:05:12.4637894Z # The new scaling point is
2022-08-20T10:05:12.4638110Z #
2022-08-20T10:05:12.4638578Z # wk := betak * sqrt(a/b) * (2*v[k]*v[k]' - J) * q
2022-08-20T10:05:12.4638818Z #
2022-08-20T10:05:12.4639082Z # with betak = W['beta'][k].
2022-08-20T10:05:12.4639302Z #
2022-08-20T10:05:12.4639636Z # 3. The scaled variable:
2022-08-20T10:05:12.4639831Z #
2022-08-20T10:05:12.4640044Z # lambda_k0 = sqrt(a*b) * c
2022-08-20T10:05:12.4640423Z # lambda_k1 = sqrt(a*b) * ( (2vk*vk' - J) * (-d*q + u/2) )_1
2022-08-20T10:05:12.4640655Z #
2022-08-20T10:05:12.4640839Z # where
2022-08-20T10:05:12.4641026Z #
2022-08-20T10:05:12.4641267Z # u = st/a - J*zt/b
2022-08-20T10:05:12.4641620Z # d = ( vk0 * (vk'*u) + u0/2 ) / (2*vk0 *(vk'*q) - q0 + 1).
2022-08-20T10:05:12.4641864Z #
2022-08-20T10:05:12.4642182Z # 4. Update scaling
2022-08-20T10:05:12.4642379Z #
2022-08-20T10:05:12.4642577Z # v[k] := wk^1/2
2022-08-20T10:05:12.4642824Z # = 1 / sqrt(2*(wk0 + 1)) * (wk + e).
2022-08-20T10:05:12.4643060Z # beta[k] *= sqrt(a/b)
2022-08-20T10:05:12.4643271Z
2022-08-20T10:05:12.4643455Z ind = m
2022-08-20T10:05:12.4643731Z for k in range(len(W['v'])):
2022-08-20T10:05:12.4643952Z
2022-08-20T10:05:12.4644314Z v = W['v'][k]
2022-08-20T10:05:12.4644524Z m = len(v)
2022-08-20T10:05:12.4644728Z
2022-08-20T10:05:12.4645040Z # ln = sqrt( lambda_k' * J * lambda_k )
2022-08-20T10:05:12.4645428Z ln = jnrm2(lmbda, n = m, offset = ind)
2022-08-20T10:05:12.4645658Z
2022-08-20T10:05:12.4645978Z # a = sqrt( sk' * J * sk ) = sqrt( st' * J * st )
2022-08-20T10:05:12.4646340Z # s := s / a = st / a
2022-08-20T10:05:12.4646607Z aa = jnrm2(s, offset = ind, n = m)
2022-08-20T10:05:12.4646898Z blas.scal(1.0/aa, s, offset = ind, n = m)
2022-08-20T10:05:12.4647123Z
2022-08-20T10:05:12.4647455Z # b = sqrt( zk' * J * zk ) = sqrt( zt' * J * zt )
2022-08-20T10:05:12.4647818Z # z := z / a = zt / b
2022-08-20T10:05:12.4648054Z bb = jnrm2(z, offset = ind, n = m)
2022-08-20T10:05:12.4648334Z blas.scal(1.0/bb, z, offset = ind, n = m)
2022-08-20T10:05:12.4648565Z
2022-08-20T10:05:12.4648857Z # c = sqrt( ( 1 + (st'*zt) / (a*b) ) / 2 )
2022-08-20T10:05:12.4649290Z cc = math.sqrt( ( 1.0 + blas.dot(s, z, offsetx = ind, offsety =
2022-08-20T10:05:12.4649582Z ind, n = m) ) / 2.0 )
2022-08-20T10:05:12.4649796Z
2022-08-20T10:05:12.4650041Z # vs = v' * st / a
2022-08-20T10:05:12.4650311Z vs = blas.dot(v, s, offsety = ind, n = m)
2022-08-20T10:05:12.4650549Z
2022-08-20T10:05:12.4651253Z # vz = v' * J *zt / b
2022-08-20T10:05:12.4651724Z vz = jdot(v, z, offsety = ind, n = m)
2022-08-20T10:05:12.4651960Z
2022-08-20T10:05:12.4652334Z # vq = v' * q where q = (st/a + J * zt/b) / (2 * c)
2022-08-20T10:05:12.4652810Z vq = (vs + vz ) / 2.0 / cc
2022-08-20T10:05:12.4653032Z
2022-08-20T10:05:12.4653355Z # vu = v' * u where u = st/a - J * zt/b
2022-08-20T10:05:12.4653655Z vu = vs - vz
2022-08-20T10:05:12.4653865Z
2022-08-20T10:05:12.4654053Z # lambda_k0 = c
2022-08-20T10:05:12.4654285Z lmbda[ind] = cc
2022-08-20T10:05:12.4654602Z
2022-08-20T10:05:12.4654873Z # wk0 = 2 * vk0 * (vk' * q) - q0
2022-08-20T10:05:12.4655466Z wk0 = 2 * v[0] * vq - ( s[ind] + z[ind] ) / 2.0 / cc
2022-08-20T10:05:12.4655699Z
2022-08-20T10:05:12.4655989Z # d = (v[0] * (vk' * u) - u0/2) / (wk0 + 1)
2022-08-20T10:05:12.4656492Z dd = (v[0] * vu - s[ind]/2.0 + z[ind]/2.0) / (wk0 + 1.0)
2022-08-20T10:05:12.4656755Z
2022-08-20T10:05:12.4657096Z # lambda_k1 = 2 * v_k1 * vk' * (-d*q + u/2) - d*q1 + u1/2
2022-08-20T10:05:12.4657651Z blas.copy(v, lmbda, offsetx = 1, offsety = ind+1, n = m-1)
2022-08-20T10:05:12.4658095Z blas.scal(2.0 * (-dd * vq + 0.5 * vu), lmbda, offset = ind+1,
2022-08-20T10:05:12.4658422Z n = m-1)
2022-08-20T10:05:12.4658904Z blas.axpy(s, lmbda, 0.5 * (1.0 - dd/cc), offsetx = ind+1, offsety
2022-08-20T10:05:12.4659242Z = ind+1, n = m-1)
2022-08-20T10:05:12.4659529Z blas.axpy(z, lmbda, 0.5 * (1.0 + dd/cc), offsetx = ind+1, offsety
2022-08-20T10:05:12.4659850Z = ind+1, n = m-1)
2022-08-20T10:05:12.4660056Z
2022-08-20T10:05:12.4660408Z # Scale so that sqrt(lambda_k' * J * lambda_k) = sqrt(aa*bb).
2022-08-20T10:05:12.4660855Z blas.scal(math.sqrt(aa*bb), lmbda, offset = ind, n = m)
2022-08-20T10:05:12.4661118Z
2022-08-20T10:05:12.4661396Z # v := (2*v*v' - J) * q
2022-08-20T10:05:12.4661843Z # = 2 * (v'*q) * v' - (J* st/a + zt/b) / (2*c)
2022-08-20T10:05:12.4662220Z blas.scal(2.0 * vq, v)
2022-08-20T10:05:12.4662532Z v[0] -= s[ind] / 2.0 / cc
2022-08-20T10:05:12.4662939Z blas.axpy(s, v, 0.5/cc, offsetx = ind+1, offsety = 1, n = m-1)
2022-08-20T10:05:12.4663341Z blas.axpy(z, v, -0.5/cc, offsetx = ind, n = m)
2022-08-20T10:05:12.4663589Z
2022-08-20T10:05:12.4663824Z # v := v^{1/2} = 1/sqrt(2 * (v0 + 1)) * (v + e)
2022-08-20T10:05:12.4664050Z v[0] += 1.0
2022-08-20T10:05:12.4664310Z blas.scal(1.0 / math.sqrt(2.0 * v[0]), v)
2022-08-20T10:05:12.4664548Z
2022-08-20T10:05:12.4664751Z # beta[k] *= ( aa / bb )**1/2
2022-08-20T10:05:12.4665097Z W['beta'][k] *= math.sqrt( aa / bb )
2022-08-20T10:05:12.4665331Z
2022-08-20T10:05:12.4665510Z ind += m
2022-08-20T10:05:12.4665708Z
2022-08-20T10:05:12.4665895Z
2022-08-20T10:05:12.4666158Z # 's' blocks
2022-08-20T10:05:12.4666360Z #
2022-08-20T10:05:12.4666633Z # Let st, zt be the updated variables in the old scaling:
2022-08-20T10:05:12.4666886Z #
2022-08-20T10:05:12.4667184Z # st = Ls * Ls', zt = Lz * Lz'.
2022-08-20T10:05:12.4667411Z #
2022-08-20T10:05:12.4667735Z # where Ls and Lz are the 's' components of s, z.
2022-08-20T10:05:12.4667992Z #
2022-08-20T10:05:12.4668302Z # 1. SVD Lz'*Ls = Uk * lambda_k^+ * Vk'.
2022-08-20T10:05:12.4668526Z #
2022-08-20T10:05:12.4668739Z # 2. New scaling is
2022-08-20T10:05:12.4668957Z #
2022-08-20T10:05:12.4669289Z # r[k] := r[k] * Ls * Vk * diag(lambda_k^+)^{-1/2}
2022-08-20T10:05:12.4669698Z # rti[k] := r[k] * Lz * Uk * diag(lambda_k^+)^{-1/2}.
2022-08-20T10:05:12.4669951Z #
2022-08-20T10:05:12.4670124Z
2022-08-20T10:05:12.4670493Z work = matrix(0.0, (max( [0] + [r.size[0] for r in W['r']])**2, 1))
2022-08-20T10:05:12.4671033Z ind = mnl + ml + sum([ len(v) for v in W['v'] ])
2022-08-20T10:05:12.4671299Z ind2, ind3 = ind, 0
2022-08-20T10:05:12.4671594Z for k in range(len(W['r'])):
2022-08-20T10:05:12.4671924Z r, rti = W['r'][k], W['rti'][k]
2022-08-20T10:05:12.4672164Z m = r.size[0]
2022-08-20T10:05:12.4672356Z
2022-08-20T10:05:12.4672564Z # r := r*sk = r*Ls
2022-08-20T10:05:12.4672858Z blas.gemm(r, s, work, m = m, n = m, k = m, ldB = m, ldC = m,
2022-08-20T10:05:12.4673126Z offsetB = ind2)
2022-08-20T10:05:12.4673545Z blas.copy(work, r, n = m**2)
2022-08-20T10:05:12.4673885Z
2022-08-20T10:05:12.4674070Z # rti := rti*zk = rti*Lz
2022-08-20T10:05:12.4674352Z blas.gemm(rti, z, work, m = m, n = m, k = m, ldB = m, ldC = m,
2022-08-20T10:05:12.4674620Z offsetB = ind2)
2022-08-20T10:05:12.4674939Z blas.copy(work, rti, n = m**2)
2022-08-20T10:05:12.4675176Z
2022-08-20T10:05:12.4696294Z # SVD Lz'*Ls = U * lmbds^+ * V'; store U in sk and V' in zk.
2022-08-20T10:05:12.4696880Z blas.gemm(z, s, work, transA = 'T', m = m, n = m, k = m, ldA = m,
2022-08-20T10:05:12.4697207Z ldB = m, ldC = m, offsetA = ind2, offsetB = ind2)
2022-08-20T10:05:12.4697652Z lapack.gesvd(work, lmbda, jobu = 'A', jobvt = 'A', m = m, n = m,
2022-08-20T10:05:12.4698006Z ldA = m, U = s, Vt = z, ldU = m, ldVt = m, offsetS = ind,
2022-08-20T10:05:12.4698319Z offsetU = ind2, offsetVt = ind2)
2022-08-20T10:05:12.4698543Z
2022-08-20T10:05:12.4698741Z # r := r*V
2022-08-20T10:05:12.4699127Z blas.gemm(r, z, work, transB = 'T', m = m, n = m, k = m, ldB = m,
2022-08-20T10:05:12.4699410Z ldC = m, offsetB = ind2)
2022-08-20T10:05:12.4699675Z blas.copy(work, r, n = m**2)
2022-08-20T10:05:12.4700021Z
2022-08-20T10:05:12.4700204Z # rti := rti*U
2022-08-20T10:05:12.4700602Z blas.gemm(rti, s, work, n = m, m = m, k = m, ldB = m, ldC = m,
2022-08-20T10:05:12.4700867Z offsetB = ind2)
2022-08-20T10:05:12.4701098Z blas.copy(work, rti, n = m**2)
2022-08-20T10:05:12.4701318Z
2022-08-20T10:05:12.4701645Z # r := r*lambda^{-1/2}; rti := rti*lambda^{-1/2}
2022-08-20T10:05:12.4702010Z for i in range(m):
2022-08-20T10:05:12.4702263Z > a = 1.0 / math.sqrt(lmbda[ind+i])
2022-08-20T10:05:12.4702566Z E ZeroDivisionError: float division by zero
2022-08-20T10:05:12.4702754Z
2022-08-20T10:05:12.4703110Z /usr/share/miniconda3/envs/test-env/lib/python3.10/site-packages/cvxopt/misc.py:628: ZeroDivisionError
Originally posted by @bnavigator in #759 (comment)
Activity
bnavigator commentedon Aug 20, 2022
Attn @Mark-Yeatman
Mark-Yeatman commentedon Aug 20, 2022
I ran the passivity tests for this branch 10 times on my local machine (Windows), couldn't reproduce. We could just remove the test case?
bnavigator commentedon Aug 20, 2022
The CI failures are on Ubuntu with conda-forge packages.
The parameter set was included for a purpose. Is it meaningful? Necessary to test?
Mark-Yeatman commentedon Aug 20, 2022
Yeah, I was trying to set something of an expectation about the relative scaling / condition number of the A and D matrices. I'll give it a look on ubuntu tonight.
Mark-Yeatman commentedon Aug 30, 2022
Looks like this hasn't reared its head again (at least according to the "actions" tab in GitHub), and I still can't reproduce locally. Do we want to close this?
murrayrm commentedon Aug 30, 2022
Actually, this error continues to be an issue. I've had to re-run most of the recent CI runs on failed tests and this was the error that showed up. I'm not sure there is a fix, but we should at least leave this open to keep track of the fact that there is an outstanding issue.
(I have to admit to being confused about why the error shows up randomly; what is the change in underlying state that leads to an error sometimes but not others?)
bnavigator commentedon Aug 30, 2022
The only difference I can imagine is that the job lands on a different worker in the GitHub Actions backend.
xfail ill-conditioned passive test (issue python-control#761)
xfail ill-conditioned passive test (issue python-control#761)
murrayrm commentedon Nov 20, 2022
I've not seen this issue pop up in a while. We can probably close unless someone has experienced a problem.
bnavigator commentedon Nov 20, 2022
The test ist marked XFAIL and still fails occasionally.
From yesterday: https://github.com/python-control/python-control/actions/runs/3502370893/jobs/5866722552#step:6:997
Mark-Yeatman commentedon Dec 12, 2022
Sorry to let this sit so long, made a pull request to fix. #814
My expectation is that now if we get a weird worker is that one of the edge cases will fail to be correctly categorized as passive instead of getting that division by zero error. ( The underlying reason could be 32bit vs 64 bit?)