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 = True2022-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 ispassive2022-08-20T10:05:12.4615096Z     return solve_passivity_LMI(sys, rho=ofp_index, nu=ifp_index) is not None2022-08-20T10:05:12.4615456Z control/passivity.py:177: in solve_passivity_LMI2022-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 sdp2022-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 conelp2022-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 variable2022-08-20T10:05:12.4621042Z         lmbda so that on exit2022-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 s2022-08-20T10:05:12.4622629Z         and z contain W^{-T}*st and W*zt, i.e, the new iterates in the current2022-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 of2022-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' blocks2022-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 = 02022-08-20T10:05:12.4627713Z         ml = len(W['d'])2022-08-20T10:05:12.4627934Z         m = mnl + ml2022-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 .* z2022-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'][:] ** -12022-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'][:] ** -12022-08-20T10:05:12.4632391Z     2022-08-20T10:05:12.4632808Z         # lmbda := s .* z2022-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_k2022-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' - J2022-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 is2022-08-20T10:05:12.4638110Z         #2022-08-20T10:05:12.4638578Z         #        wk := betak * sqrt(a/b) * (2*v[k]*v[k]' - J) * q2022-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) * c2022-08-20T10:05:12.4640423Z         #        lambda_k1 = sqrt(a*b) * ( (2vk*vk' - J) * (-d*q + u/2) )_12022-08-20T10:05:12.4640655Z         #2022-08-20T10:05:12.4640839Z         #    where2022-08-20T10:05:12.4641026Z         #2022-08-20T10:05:12.4641267Z         #        u = st/a - J*zt/b2022-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 scaling2022-08-20T10:05:12.4642379Z         #2022-08-20T10:05:12.4642577Z         #        v[k] := wk^1/22022-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 = m2022-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 / a2022-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 / b2022-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 / a2022-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 / b2022-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 / cc2022-08-20T10:05:12.4653032Z     2022-08-20T10:05:12.4653355Z             # vu = v' * u  where u =  st/a - J * zt/b2022-08-20T10:05:12.4653655Z             vu = vs - vz2022-08-20T10:05:12.4653865Z     2022-08-20T10:05:12.4654053Z             # lambda_k0 = c2022-08-20T10:05:12.4654285Z             lmbda[ind] = cc2022-08-20T10:05:12.4654602Z     2022-08-20T10:05:12.4654873Z             # wk0 = 2 * vk0 * (vk' * q) - q02022-08-20T10:05:12.4655466Z             wk0 = 2 * v[0] * vq - ( s[ind] + z[ind] ) / 2.0 / cc2022-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/22022-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, offsety2022-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, offsety2022-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) * q2022-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 / cc2022-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.02022-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/22022-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 += m2022-08-20T10:05:12.4665708Z     2022-08-20T10:05:12.4665895Z     2022-08-20T10:05:12.4666158Z         # 's' blocks2022-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 is2022-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, 02022-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*Ls2022-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*Lz2022-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*V2022-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*U2022-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 zero2022-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