Skip to content

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

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

Description

@Mar-tian
Mar-tian
opened on 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));

Activity

  1. murrayrm commented on Jan 7, 2020

    @murrayrm
    murrayrm
    on Jan 7, 2020
    Member

    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?

  2. roryyorke commented on Jan 7, 2020

    @roryyorke
    roryyorke
    on Jan 7, 2020
    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.

  3. Mar-tian commented on Apr 2, 2020

    @Mar-tian
    Mar-tian
    on Apr 2, 2020
    Author

    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?

  4. repagh commented on May 5, 2020

    @repagh
    repagh
    on May 5, 2020
    Member

    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?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

      Milestone

      No milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions