Skip to content

Removed epsilon perturbation value in solve_passivity_LMI. Fix associated unit test. #814

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 3 commits into from
Dec 17, 2022
Merged

Removed epsilon perturbation value in solve_passivity_LMI. Fix associated unit test. #814

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
4792163
Removed epsilon perturbation value in solve_passivity_LMI. Change uni…
Mark-Yeatman Dec 12, 2022
a17e3c8
Remove test case.
Mark-Yeatman Dec 15, 2022
e6b837d
Update control/passivity.py
Mark-Yeatman Dec 16, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
31 changes: 17 additions & 14 deletions control/passivity.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,6 @@
except ImportError:
cvx = None

eps = np.nextafter(0, 1)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

(This is a by-the-way comment, feel free to ignore.)

This eps was very small; I think a more typical "eps" would be "the smallest number eps such that 1+eps > 1", or a multiple or maybe square root of that. np.finfo is useful to find these sorts of value.

In [402]: np.nextafter(0, 1)
Out[402]: 5e-324

In [403]: np.finfo(float).eps
Out[403]: 2.220446049250313e-16

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That is a pro tip, and what I was aiming for.


__all__ = ["get_output_fb_index", "get_input_ff_index", "ispassive",
"solve_passivity_LMI"]
Expand All @@ -28,8 +27,8 @@ def solve_passivity_LMI(sys, rho=None, nu=None):
passive. Inputs of None for `rho` or `nu` indicate that the function should
solve for that index (they are mutually exclusive, they can't both be
None, otherwise you're trying to solve a nonconvex bilinear matrix
inequality.) The last element of `solution` is either the output or input
passivity index, for `rho` = None and `nu` = None respectively.
inequality.) The last element of the output `solution` is either the output or input
passivity index, for `rho` = None and `nu` = None respectively.

The sources for the algorithm are:

Expand Down Expand Up @@ -74,11 +73,8 @@ def solve_passivity_LMI(sys, rho=None, nu=None):
D = sys.D

# account for strictly proper systems
[n, m] = D.shape
D = D + eps * np.eye(n, m)

[_, m] = D.shape
[n, _] = A.shape
A = A - eps*np.eye(n)

def make_LMI_matrix(P, rho, nu, one):
q = sys.noutputs
Expand Down Expand Up @@ -113,7 +109,7 @@ def make_P_basis_matrices(n, rho, nu):
P[i, j] = 1
P[j, i] = 1
matrix_list.append(make_LMI_matrix(P, 0, 0, 0).flatten())
zeros = eps*np.eye(n)
zeros = 0.0*np.eye(n)
if rho is None:
matrix_list.append(make_LMI_matrix(zeros, 1, 0, 0).flatten())
elif nu is None:
Expand Down Expand Up @@ -149,9 +145,9 @@ def P_pos_def_constraint(n):
if rho is not None and nu is not None:
sys_constants = -make_LMI_matrix(np.zeros_like(A), rho, nu, 1.0)
elif rho is not None:
sys_constants = -make_LMI_matrix(np.zeros_like(A), rho, eps, 1.0)
sys_constants = -make_LMI_matrix(np.zeros_like(A), rho, 0.0, 1.0)
elif nu is not None:
sys_constants = -make_LMI_matrix(np.zeros_like(A), eps, nu, 1.0)
sys_constants = -make_LMI_matrix(np.zeros_like(A), 0.0, nu, 1.0)

sys_coefficents = np.vstack(sys_matrix_list).T

Expand All @@ -174,8 +170,15 @@ def P_pos_def_constraint(n):

# crunch feasibility solution
cvx.solvers.options['show_progress'] = False
sol = cvx.solvers.sdp(c, Gs=Gs, hs=hs)
return sol["x"]
try:
sol = cvx.solvers.sdp(c, Gs=Gs, hs=hs)
return sol["x"]

except ZeroDivisionError as e:
raise ValueError("The system is probably ill conditioned. "
"Consider perturbing the system matrices by a small amount."
) from e



def get_output_fb_index(sys):
Expand All @@ -194,7 +197,7 @@ def get_output_fb_index(sys):
float
The OFP index
"""
sol = solve_passivity_LMI(sys, nu=eps)
sol = solve_passivity_LMI(sys, nu=0.0)
if sol is None:
raise RuntimeError("LMI passivity problem is infeasible")
else:
Expand All @@ -218,7 +221,7 @@ def get_input_ff_index(sys):
float
The IFP index
"""
sol = solve_passivity_LMI(sys, rho=eps)
sol = solve_passivity_LMI(sys, rho=0.0)
if sol is None:
raise RuntimeError("LMI passivity problem is infeasible")
else:
Expand Down
10 changes: 4 additions & 6 deletions control/tests/passivity_test.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -99,16 +99,14 @@ def test_system_dimension():


@pytest.mark.parametrize(
"systemmatrices, expected",
[((A, B, C, D*0.0), True),
"system_matrices, expected",
[((A, B, C, D*0), True),
((A_d, B, C, D), True),
pytest.param((A*1e12, B, C, D*0), True,
marks=pytest.mark.xfail(reason="gh-761")),
((A, B*0, C*0, D), True),
((A*0, B, C, D), True),
((A*0, B*0, C*0, D*0), True)])
def test_ispassive_edge_cases(systemmatrices, expected):
sys = ss(*systemmatrices)
def test_ispassive_edge_cases(system_matrices, expected):
sys = ss(*system_matrices)
assert passivity.ispassive(sys) == expected


Expand Down