misconstructed plant causes dkpy to hang indefinitely using 100% cpu for hours with no progress

[bbbccd]

Expected Behavior

dkpy should detect malformed plants if practical, or crash with helpful messages, or provide some way to see (lack of) progress. (I will open a separate issue for the latter). ctrl+c should kill the program.

Actual Behavior

A plant was created missing one of the middle outputs. This could have been identified by a mismatch between the plant output count and (ny+blocksize). Instead it caused dkpy to hang indefinitely using 100% cpu for hours with no progress. it did not respond to ctrl+c. had to kill with killall -9 python

Steps to Reproduce the Problem

Let me know if you want the code with the bad plant. it's ~200 lines of code. I haven't tried to create a minimal example, but I can work on that if you prefer that over seeing the original code.

Specifications

• Package version: 0.1.8
• Python version: Python 3.13.3
• Platform: arch

[sdahdah]

Thanks for reporting. Please post the code and we will take a look

[bbbccd]

what's the best way to post long code? I have just put it here in a giant code block. should i attach a file?

from control import (ss, zpk, lqr, lqe, c2d, sample_system, hinfsyn, series,
                    tf, bode_plot, hankel_singular_values, balanced_reduction,
                    singular_values_plot)
from numpy import block, diag, array, zeros, logspace
from math import sqrt
from matplotlib import pyplot as plt
import dkpy
#
# model.m was converted into model.py and abandoned
# model.py is thus more recent

# also of note:
# length = 0.553m
# center of axle to bottom of board = 0.047m
#
# when free rolling horizontal, my food scale reads this weight
# 
#    pwm    weight
#   -255    250g
#   -200    150g
#   -100    150g
#    -50    15g
#     50    10g
#    100    200g
#    200    200g
#    255    260g
# 
# which is incredibly nonlinear with a saturation torque of 0.184N-m
# gear jamming contributes about 100-200g in either direction :o :)
# that's like 100% full scale power
# thats so awful :)
#
# free spinning (airborn) wheels spin on average 1.2m/s forwards and 1.0m/s backwards
#
# 8 tooth motor gear meshes with 74 tooth gear
# the 74 tooth gear is bonded to a 15 tooth gear that meshes with a 47 tooth gear
# and so on. so the gear ratio is (74/8)*(47/15)*(47/25)=54.489
#    8*
#     #                                                                         
#     #=====15                                                                  
#  ---#=====--------                                                            
#     #=====                                                                    
#   74#  #                                                                      
#        #=====25                                                             
#        #=====                                                                 
#  ------#=====-----                                                            
#        #=====                                                                 
#        #=====                                                                 
#      47#  #                                                                   
#           #=====25                                                            
#           #=====                                                              
#  ---------#=====--                                                            
#           #=====                                                              
#           #=====                                                              
#         47#                                                                   
#                                                                               
#                                                                               

r=0.075                 # [m]       wheel radius
m =  0.394              # [kg]      wheel and axle mass
M = 1.574               # [kg]      everything else with all the stuff on it
I_w = 0.5*0.218*r**2    # [kg-m**2]  0.218 is mass of both wheels without axle
L = 0.19                # [m]       distance from axle to center of mass of the board with all the stuff on it
p = 4.6                 #           experimentally observed natural pole of theta/tau
g = 9.6                 # [m/s**2]   gravity



# the experimentally measured natural pole of theta/tau is 4.6
# solving for the moment of inertia of the board then gives 0.1219kg-m**2
# if we assume the center of mass is in the middle of the board...
# (it's not because of all the stuff on it)
# ...and evenly distributed (again, no) we get 0.0397kg-m**2, a third
# i use the geometric average of these values
# which gives a natural pole of 5.9
I1=L*M*(-L*I_w*p**2 - L*m*p**2*r**2 + M*g*r**2 + g*I_w + g*m*r**2)/(p**2*(M*r**2 + I_w + m*r**2))
I2=1/12*M*.55**2
I=sqrt(I1*I2)
p = sqrt((L*M**2*g*r**2 + L*M*g*I_w + L*M*g*m*r**2)/(I*M*r**2 + I*I_w + I*m*r**2 + L**2*M*I_w + L**2*M*m*r**2))



# state space model
# x = [theta thetad x xd]
# y = [theta x]'
A = array([ [0                                                                                          ,1  ,0  ,0],
            [L*M*g*(M*r**2 + I_w + m*r**2)/(I*M*r**2 + I*I_w + I*m*r**2 + L**2*M*I_w + L**2*M*m*r**2)   ,0  ,0  ,0],
            [0                                                                                          ,0  ,0  ,1],
            [L**2*M**2*g*r**2/(I*M*r**2 + I*I_w + I*m*r**2 + L**2*M*I_w + L**2*M*m*r**2)                ,0  ,0  ,0]     ])
B = array([ [0],
            [(L*M*r + M*r**2 + I_w + m*r**2)/(I*M*r**2 + I*I_w + I*m*r**2 + L**2*M*I_w + L**2*M*m*r**2)],
            [0],
            [r*(I + L**2*M + L*M*r)/(I*M*r**2 + I*I_w + I*m*r**2 + L**2*M*I_w + L**2*M*m*r**2)]                 ])
C = array([ [1  ,0  ,0  ,0],
            [0  ,0  ,1  ,0]])
D = array([[0],[0]])

#rlocus(ss(A,B,C)(2,1))

# h infinity model
# w = [unmodeled_torque board wheel angle_noise, position_error]
# z = [theta x u]
B_disturbance = array([     [0      ,0      ,0      ,0      ,0],
                            [B[1,0] ,B[1,0] ,0      ,0      ,0],
                            [0      ,0      ,0      ,0      ,0],
                            [B[3,0] ,0      ,B[3,0] ,0      ,0]     ])
D_actuator_cost = array([[0],[0],[1]]);
D_noise = array([   [0  ,0  ,0  ,1  ,0],
                    [0  ,0  ,0  ,0  ,1]     ])
P = ss(     A,
block(  [   [B_disturbance,B]   ]),
block(  [   [C],
            [0,0,0,0],
            [C]]),
block(  [   [zeros((3,5)) ,D_actuator_cost], 
            [D_noise    ,zeros((2,1))]])
);


# weights
# 
# for input weights, i just choose a "typical" disturbance
# 
# 100g gear jam measured by scale in experiment above -> 0.98N -> 0.0735N-m
# and we'll just chose the same for board-only and wheel-only disturbances
# 5 degree angle measurement error
# backlash + one encoder count position error = 13mm+2mm = 15mm
#
# for output weights, i choose "equal terribleness"
# 10 degrees, 50cm, control saturation 
#
# y and u pass through unweighted
W_input = diag([0.0735 ,0.0735 ,0.0735 ,0.087266 ,0.015 ,1]);
W_output = ss(tf([[[1/.1745] ,[0]    ,[0]               ,[0]    ,[0]    ], 
                [    [0]     ,[1/0.5],[0]               ,[0]    ,[0]    ], 
                [    [0]     ,[0] ,[1e6/.184,5e6/.184]  ,[0]    ,[0]    ], 
                [    [0]     ,[0]    ,[0]               ,[1]    ,[0]    ], 
                [    [0]     ,[0]    ,[0]               ,[0]    ,[1]    ]], 
               [[    [1]     ,[1]    ,[1]               ,[1]    ,[1]    ], 
                [    [1]     ,[1]    ,[1]               ,[1]    ,[1]    ], 
                [    [1]     ,[1]  ,[1,5e6]             ,[1]    ,[1]    ], 
                [    [1]     ,[1]    ,[1]               ,[1]    ,[1]    ], 
                [    [1]     ,[1]    ,[1]               ,[1]    ,[1]    ]]))
P = W_output * P * W_input;




K_hinf, N, gamma, info = hinfsyn(P,2,1);
#K = balanced_reduction(K,4)
#bode_plot([K[0,0],balanced_reduction(K,4)[0,0]])
#singular_values_plot([P.lft(K),P.lft(balanced_reduction(K,4))])
#plt.show()
K_hinf = balanced_reduction(K_hinf,4)
#bode_plot(K_hinf[0,0])
#plt.show()

# LQG controller
#
# we use bryson's rule for lqr weights, using similar worst case scenarios as above
# we'll make worst case theta_d to be worst case theta in a quarter second
# and worst case x_d to be half full speed
# 
# we'll invert that rule :) for lqe weights. don't know dual to brysons rule...
# theta and x are only disturbed through their derivative states,
# so we make their variances small
# 
sysss = ss(A,B,C,D);
K_lqr,_,_=lqr(sysss,diag([1/0.1745**2 ,1/(0.1745*4)**2 ,1/0.5**2 ,1/0.5**2]),diag([1/0.184**2]));
#print(diag([0.0001**2 ,(B[1,0]*0.0735)**2 ,0.0001**2 ,(B[3,0]*0.0735)**2]).shape)
#L=lqe(sysss,diag([0.0001**2 ,(B[1,0]*0.0735)**2 ,0.0001**2 ,(B[3,0]*0.0735)**2]),diag([0.087266**2,0.015**2]));
L,_,_=lqe(A,diag([1,1,1,1]),C,diag([0.0001**2 ,(B[1,0]*0.0735)**2 ,0.0001**2 ,(B[3,0]*0.0735)**2]),diag([0.087266**2,0.015**2]));

Clqg=ss(sysss.A-sysss.B@K_lqr-L@sysss.C,L,K_lqr,[[0,0]]);
Nlqg=P.lft(-Clqg,1,2);


#singular_values_plot([P.lft(K_hinf),P.lft(Clqg)])
#singular_values_plot(P.lft(Clqg))
#singular_values_plot(P.lft(K_hinf))
#plt.show()

# THE RESULT!!!! :)
K_arduino = sample_system(K_hinf,0.0025);

# Mu Synthesis / DK Iteration
# outputs: theta theta theta x u(weighted) theta x
#                     |<----z------------>|<--y-->|
# 
# inputs:  theta*delta theta*delta unmodeled_torque board wheel angle_noise position_error u
#                                 |<------------------------w---------------------------->|u
# states:  theta thetad x xd
P = ss(     A,
block(  [[      array([[0],[0.2*A[1,0]],[0],[0]])   ,array([[0],[0],[0],[0.2*A[3,0]]])   ,B_disturbance  ,B       ]]),
block(  [   [1,0,0,0],                           # theta
            [1,0,0,0],                           # theta
            [C],                                # theta x
            [0,0,0,0],                          # u
            [C]]),                              # theta x
block(  [   [zeros((2,8))                                  ],
            [zeros((3,2))  ,zeros((3,5))   ,D_actuator_cost], 
            [zeros((2,2))  ,D_noise        ,zeros((2,1))   ]])
);

dk_iter = dkpy.DkIterListOrder(
    controller_synthesis=dkpy.HinfSynSlicot(),
    structured_singular_value=dkpy.SsvLmiBisection(),
    d_scale_fit=dkpy.DScaleFitSlicot(),
    fit_orders=[4, 4, 4],)
block_structure = array([[-1, 0], [-1, 0], [5, 3]]) # NotImplementedError: Real perturbations are not yet supported.
block_structure = array([[1, 0], [1, 0], [5, 3]]) # NotImplementedError: Diagonal perturbations are not yet supported.
block_structure = array([[1, 1], [1, 1], [5, 3]]) # NotImplementedError: Nonsquare perturbations are not yet supported.
P = array([[  1   ,0  ,0  ,0  ,0  ,0  ,0  ],
     [  0   ,1  ,0  ,0  ,0  ,0  ,0  ],
     [  0   ,0  ,1  ,0  ,0  ,0  ,0  ],
     [  0   ,0  ,0  ,1  ,0  ,0  ,0  ],
     [  0   ,0  ,0  ,0  ,0  ,0  ,0  ],
     [  0   ,0  ,0  ,0  ,0  ,0  ,0  ],
     [  0   ,0  ,0  ,0  ,0  ,1  ,0  ],
     [  0   ,0  ,0  ,0  ,0  ,0  ,1  ]])*P # add ficticious identically zero performance outputs to play well with square performance uncertainty block
block_structure = array([[1, 1], [1, 1], [5, 5]])
K, N, mu, d_scale_fit_info, info = dk_iter.synthesize(
    P,2,1,logspace(-3, 3, 61),
    block_structure,
)

[sdahdah]

Code block is fine, thanks. You can add syntax highlighting using

```python

as the starting line