Repository navigation
hinfsyn in python-control is returning different results as compared to hinfsyn in MATLAB's robust control toolbox #367
Description
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
murrayrm commentedon Jan 7, 2020
I tried this example out in Python 3.7, latest master and the call to
hinfsyndidn't return after running for 10 minutes (at which point I lost patience). Thehinfsynfunction is a pass through toslycotusing thesb10adfunction.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 functiontd04ad.As a check, I entered the state space realization that MATLAB compute. No difference.
I next tried calling
slycot.sb10addirectly using the matrices from the example above:This is the same call as
hinfsynmakes, except that I used an initial guess ofgamma=10instead ofgamma=1e100(value insidehinfsyn, which is probably not the best guess to use). This returned very quickly with the errorBottom line: nothing wrong with your Python code, but something wrong in
slycot/sb10ad.@roryyorke @repagh @bnavigator Any clues?
roryyorke commentedon Jan 7, 2020
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 commentedon Apr 2, 2020
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 commentedon 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 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?