Skip to content

hinfsyn in python-control is returning different results as compared to hinfsyn in MATLAB's robust control toolbox #367

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

Jump to bottom
Open
Mar-tian opened this issue Jan 7, 2020 · 4 comments
Open

hinfsyn in python-control is returning different results as compared to hinfsyn in MATLAB's robust control toolbox #367

Mar-tian opened this issue Jan 7, 2020 · 4 comments

Comments

@Mar-tian
Copy link

Mar-tian commented Jan 7, 2020

Hi,

I am very new to control systems, and coding in general. I am currently exploring the H-infinity robust control method and tried implementing the hinfsyn function on a simple spring mass damper system.

The following is the code I'm trying to implement:

import` control
import numpy


# spring mass damper system
m = 1  # 1 kg
b = 12 # 12 N-s/m
k = 20 # 20 N/m


num = [[[-1], [-1]],[[0], [1]], [[-1], [-1]]]
den = [[[1,1], [m, b, k]], [[1], [1]], [[1,1], [m, b, k]]]
Plant = control.tf2ss(num, den)
#print(Plant)


nmeas = 1
ncont = 1
K, CL, gamma, rcond = control.hinfsyn(Plant, nmeas, ncont)
print(K)
print(gamma)
Ktf = control.ss2tf(K);
print(Ktf);`

The code returns "A stabilizing controller cannot be found"

However, a similar code returns a controller on MATLAB. I have checked and the state space variable "Plant" returned by python-control.tf2ss and MATLAB's tf2ss is almost the same (I'm not sure why such a difference arises but I'm assuming MATLAB scales the matrices A, B, C, D).

Could someone help me out, I'm not sure what I'm doing wrong with my python code. Also let me know if this the wrong forum for posting such a query. I wasn't sure where else to go.

For reference:

I was implementing the model described in Robust Control Design with MATLAB by D.-W. Gu, P. Hr. Petkov and M. M. Konstantinov (Pages 36, 37)

The following was my corresponding code in MATLAB:

% spring mass damper system
m = 1;  % 1 kg
b = 12; % 12 N-s/m
k = 20; % 20 N/m

% plant transfer function 
Gnom = [-tf(1,[1,1]) -tf(1,[m, b, k]);tf(0,1) tf(1,1); -tf(1,[1,1]) -tf(1,[m, b, k])]; 

% state space model of the plant method
P=ss(Gnom); 

% applying hinfinity
nmeas = 1;
ncont = 1;
[K, CL, gamma] = hinfsyn(P, nmeas, ncont);
disp(gamma);

% displaying the state space model properties of the controller
disp('K is');
disp(K);
disp(K.A);
disp(K.B);
disp(K.C);
disp(K.D);
disp(size(K));

%retrieving transfer function of controller method

[numk, denk] = ss2tf(K.A, K.B, K.C, K.D); 
Ktf = tf(numk, denk);
disp('Ktf is');
disp(Ktf);
disp(size(Ktf));
@murrayrm
Copy link
Member

murrayrm commented Jan 7, 2020

I tried this example out in Python 3.7, latest master and the call to hinfsyn didn't return after running for 10 minutes (at which point I lost patience). The hinfsyn function is a pass through to slycot using the sb10ad function.

The issue with the state space dynamics being different in python-control versus MATLAB likely has to do with the way that python-control does transformation to state space. This is also a pass through to slycot, using the function td04ad.

As a check, I entered the state space realization that MATLAB compute. No difference.

I next tried calling slycot.sb10ad directly using the matrices from the example above:

A = [[-1, 0 , 0], [0, -12, -5], [0, 4, 0]] 
B = [[2,0], [0, 0.5], [0, 0]]   
C = [[-0.5, 0, -0.5], [0, 0, 0], [-0.5, 0, -0.5]] 
D = [[0, 0],  [0, 1], [0, 0]]
out = slycot.sb10ad(3, 2, 3, 1, 1, 10, A, B, C, D)

This is the same call as hinfsyn makes, except that I used an initial guess of gamma=10 instead of gamma=1e100 (value inside hinfsyn, which is probably not the best guess to use). This returned very quickly with the error

ValueError: A stabilizing controller cannot be found.

Bottom line: nothing wrong with your Python code, but something wrong in slycot/sb10ad.

@roryyorke @repagh @bnavigator Any clues?

@roryyorke
Copy link
Contributor

No, sorry. I imagine Matlab's H-infinity solver is better than our oldish (I think?) SLICOT version; MathWorks probably license the latest SLICOT, or have other algorithms altogether.

It's a bit odd that @Mar-tian doesn't get the same "[h-]infinite wait" (sorry) a la #152 that you did.

@Mar-tian
Copy link
Author

Mar-tian commented Apr 2, 2020 •
edited
Loading

Hello,

@murrayrm and @roryyorke. Thank you for your reply. I somehow never received a notification that you guys replied to the issue back then and would like to apologise for not getting back. Thank you once again.

I have only recently got back to trying to solve this problem on Python again and realized that you guys had replied back.

Anyway, over the past few months, I went back to working on this problem on MATLAB and while I might not be able to contribute much to the solution, I would like to update you guys.

At the time of posting this problem, even though the hinfsyn function in MATLAB returns a gamma value (something in the order of 10e-5) and some controller, the controller is basically not usable as the numerator of the controller transfer function turns out be 0. What solved this problem was using the augw function (python's control toolbox also has a similar function) which enables me to add certain weighting functions to my "P" state space system above (shall refrain from going into further details as they are not really relevant). This returned a gamma value of approximately 1 and a usable controller state space model/ transfer function.

Further, MATLAB also has a way of manipulating the gamma range in the arguments of the hinfsyn function, which also returns a usable controller.

Today, while trying to move my MATLAB code to python, I used the augw function in python's control toolbox in my code, and still my python code ran for a long time (probably half an hour?) and returned "stabilizing controller cannot be found". I shall try changing the initial gamma value in slycot.sb10ad and see if it makes any difference.

I'm not really familiar with the control toolbox in python yet, but would you be knowing of anyway to limit my gamma to some particular range (say 10 - 1000)?

Thank you once again for replying.

P.S. I hope you guys are healthy and safe in these times.

Update:

Have tried both the suggestions to my new code (changing the initial value of gamma in robust.py and calling slycot.sb10ad()). Both of them are still returning "stabilizing controller cannot be found" but in a tolerable time.

Btw, could you clarify something. Setting the gamma value in robust.py sets the upper bound of gamma or the lower bound?

@repagh
Copy link
Member

repagh commented May 5, 2020

I tried another one of the Slycot routines, sb10fd. I have that currently in a Slycot branch, which I can create a pull request for when we have the new exception framework in place.

It also fails, with a specific error; number 8,

                The Y-Riccati equation was not solved
                successfully (the controller is not admissible or
                there are numerical difficulties)

The system also looks a bit fishy, with the second output having no dynamics, just a direct feedthrough from the second input. Are you sure this is the correct system?

@roryyorke roryyorke mentioned this issue May 17, 2020
Closed
@NaoPross NaoPross mentioned this issue Feb 29, 2024
Closed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
4 participants
@roryyorke @murrayrm @repagh @Mar-tian