Passivity test is ill-conditioned
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