I get a robust controller result using musyn that is empty; "Krob = [ ]".

Below is a Live Editor script of a simple mass/spring/damper system that I'm trying to design a robust controller using the musyn command. This work follows the example in Zhou and Doyle, Essentials of Robust Control, and www.youtube.com/@artcellCTRL. I am unable to figure out why I get the "Krob = []" output (see the output at the end of the script). The 'K Step', 'Peak MU', 'D Fit' are Inf, and the 'D' result is NaN in the D-K interation. I don't understand what the problem is.

Thanks so much.

clear all;close all;clc;

Simple spring/mass/damper system is described as follows:

This system can be represented by the following differential equation of motion:

and represented by the following continuous time block diagram and state-space representation:

with m = mass of the block, k = spring constant and c = damping coefficient.

We wish to represent the parameters m, k and c with uncertainty using Eta and Delta. Where Delta will be set to vary between -1 to +1, and Eta will set the limit of the total variation (e.g., maximum percent variation) using the 'ureal' command:

delta_m = ureal("delta_m",1,'PlusMinus',1);% Uncertain LTI dynamics, +/-1
delta_k = ureal("deltal_k",1,'PlusMinus',1);% Uncertain LTI dynamics, +/-1 
delta_c = ureal("delta_c",1,'PlusMinus',1);% Uncertain LTI dynamics, +/-1 
%

Setting Eta at 10% maximum variation:

eta_m=0.1;% vary 10%
eta_k=0.1;% vary 10%
eta_c=0.1;% vary 10%

Re-drawing the block diagram to include the variation of the parameters:

Expanding the block diagram for c and k:

Next, we label all of the block connections along with defining the dummy block ap1 to ensure each block I/O have a unique identifier and labeling the 1/s blocks a1 and a2:

ap1=ss(1);
ap1.u='s1';
ap1.y='yx1';

a1=tf([1],[1 0]);% 1/s term
a1.u='s2';
a1.y='s1';
%
a2=tf([1],[1 0]);% 1/s term
a2.u='s8';
a2.y='s2';

Sum1=sumblk('s7=u-s9');%CORRECTED PER RESPONSE TO POST
Sum2=sumblk('s9=s4+s3');

Next, we define the LFT of the uncertain parameter blocks:

^^^CORRECTED PER RESPONSE TO POST^^^

^^^CORRECTED PER RESPONSE TO POST^^^

(See Zhou and Doyle, Essentials of Robust Control pp166)

Next, we define each LFT block and code into Matlab:

MASS:

m0=1;
Mm=[-eta_m, 1;-eta_m/m0, 1/m0 ]; % 
Mm_sys=ss(Mm);
Mm_sys.u={'um','s7'};
Mm_sys.y={'ym','s8'};

SPRING CONSTANT, k:

^^^CORRECTED PER RESPONSE TO POST^^^

k0=1;
%Mk=[k0*eta_k, 0;k0, 1];% ORIGINAL
Mk=[0, k0*eta_k;1, k0];%CORRECTED PER RESPONSE TO POST
Mk_sys=ss(Mk);
Mk_sys.u={'uk','s1'};
Mk_sys.y={'yk','s3'};

DAMPING COEFFICIENT, c0:

^^^CORRECTED PER RESPONSE TO POST^^^

c0=1;
%Mc=[c0*eta_c, 0;c0, 1];% ORIGINAL
Mc=[0, c0*eta_c;1, c0];%CORRECTED PER RESPONSE TO POST
Mc_sys=ss(Mc);
Mc_sys.u={'uc','s2'};
Mc_sys.y={'yc','s4'};

Next, we re-draw the LFT block diagram and set up the operation of 'Pulling Out the Deltas' for the Mu Synthesis architecture:

Pulling Out the Deltas:

Pulling out the deltas, forming the delta matrix block and defining the I/O for the mu synthesis architecture:

Adding a reference input, r and generating an error output of the position state, x1, along with pulling out the velocity state x2 of the system.

%
Sum3=sumblk('ye=r-yx1');
%

ap2=ss(1);
ap2.u='s2';
ap2.y='yx2';

ap3=ss(1);
ap3.u='yx1';
ap3.y='yx1_o';

Now we form the delta block:

delta_mtx=[delta_m, 0, 0;0, delta_c, 0;0 , 0, delta_k];
delta_sys=ss(delta_mtx);
delta_sys.u={'ym','yc','yk'};
delta_sys.y={'um','uc','uk'};

Next, we connect all of the components to form the Mu Synthesis architecture, adding, for reference, the feedback K block that will be calculated with the musyn command later on:

P1=connect(a1,a2,Mm_sys,Mk_sys,Mc_sys,Sum1,Sum2,ap1,ap2,delta_sys,Sum3, {'um','uc','uk','u','r'},{'ym','yc','yk','ye','u','ye','yx2'});

For reference, the input/output of the Mu Synthesis block in the Matlab 'connect' command is as follows:

Next, to simplify, we ignore the delta terms and allow them to be lumped into the final Mu Synthesis block so we only have w = {'r'}, and u = {'u'} as inputs and z = {'ye', 'u'}, and y = {'ye', 'yx2'} as the outputs to form the final Mu Synthesis block P2:

Now we re-write the connect command for P2 (which is the same as the one for P1, but with w = {'r'}, and u = {'u'} as inputs and z = {'ye', 'u'}, and y = {'ye', 'yx2'} as the outputs):

%P2=connect(a1,a2,Mm_sys,Mk_sys,Mc_sys,Sum1,Sum2,ap1,ap2,delta_sys,Sum3,{'u','r'},{'ye','u','ye','yx2'});% ORIGINAL
P2=connect(a1,a2,Mm_sys,Mk_sys,Mc_sys,Sum1,Sum2,ap1,ap2,delta_sys,Sum3,{'r','u'},{'ye','u','ye','yx2'});%CORRECTED PER RESPONSE TO POST (THIS WAS THE KEY TO MY ERROR!!)

The resultant connect command I/O and block function is as follows:

CORRECTED PER RESPONSE TO QUESTION (swapped 'u' and 'r')

V V V V V

Finally, we setup the musyn command by defining the number of outputs, 'ncont' and the number of measured states as inputs, 'nmeas' of the controller, Krob and executing the musyn command:

%
ncont = 1; % one control signal, 'u'
nmeas = 2; % 2 measurement signals,'ye' and 'yx2'
%

[Krob,CLperf] = musyn(P2,nmeas,ncont)

ORIGINAL POST QUESTION:

• After running the script, I get the output shown below. I don't understand why. I tried changing the mass, damper and spring constant parameters, but I still get the same output.
• Any ideas as to what is going on?

************

D-K ITERATION SUMMARY:

-----------------------------------------------------------------

Robust performance Fit order

-----------------------------------------------------------------

Iter K Step Peak MU D Fit D

1 Inf Inf Inf NaN

Best achieved robust performance: Inf

Krob =

[]

************

************

************

THE CORRECTIONS PROVIDED IN RESPONSE TO THE POST FIXED THE PROBLEM!

NEXT STEPS:

1) Figuring out how to define and apply the weights.

2) Interpreting and figuring out how to implement the results.

[Paul]

Hi Peter,

I see that you've edited the question again, but I have no idea what the edit was or if further response is needed. If the edit(s) is something substantive, i.e., more than just correcting a typo or clarifying the wording, it woud be immensely helpful to you to add a comment in the answer thread below that explains the update and provides insight into any current issues, should there be any. Good luck with the rest of your project if you've got it all sorted out.

Answer 1

[Paul]

Hi Peter,

My compliments on the quality of and effort put into this post!

Based on only a quick inspection, I have the following comments:

1. c0 = k0 = 1, but the sum1 and sum2 are implemented as positive summations. The system will therefore be unstable, which isn't necessarily a problem, but seems like might not be what's intended to be modeled for a basic mass-spring-damper system.

2. The LFT representations for the uncertain spring and damper appear to be incorrect. I did not check the LFT representation for the mass. For the dashpot, we can see by inspection

s4 = c0*s2 + uc;

yc = etac*c0*s2

uc = deltac*yc

Or in matrix form

[yc ; s4] = [0 etac*c0 ; 1 c0 ] * [uc ; s2]

uc = deltac*yc

The graphical representations of the LFTs for k and c have the arrows going left to right; better to have them go right to left as for the LFT for mass.

3. Shouldn't the performance output z = [ye ; u] mulitplied by frequency dependent weights to drive the controller to achieve the desired shapes for the transfer functions Ye(s)/R(s) and U(s)/R(s)? Typically the weights will be selected to make the former small at low frequency and the latter small at high frequency. The function makeweight can be used to create the weights in the frequency domain; see Wl and Wh on that doc page. An example at musyn shows how to define and apply a weight to the error transfer function, aka, the output sensitivity function.

4. The text of the post suggest the uncertainty on the parameters should be 0-10%, but it looks like the model is implementing +-10%.

[Paul]

I was so enthralled by the all of the graphics and equations in the question that I missed the fact that there isn't a reason to explicitly do all of the LFT work to model the uncertain parameters. Instead, just let the toolbox do the work. Here's how the plant could be built with +-10% uncertainty on m, c, and k (and assuming the nominal plant is the standard mass-spring-damper: M*yddot + C*ydot + K*y = F )

m0 = 1;
c0 = 1;
k0 = 1;
m = ureal("m",m0,'Percentage',10);
c = ureal("c",c0,'Percentage',10);
k = ureal("k",k0,'Percentage',10);

Edit: original code was missing 1/m in the B matrix.

%plant = ss([0 1;-k/m -c/m],[0;1],eye(2),0,'InputName','u','OutputName',{'x1','x2'})
plant = ss([0 1;-k/m -c/m],[0;1/m],eye(2),0,'InputName','u','OutputName',{'x1','x2'})

Uncertain continuous-time state-space model with 2 outputs, 1 inputs, 2 states.
The model uncertainty consists of the following blocks:
  c: Uncertain real, nominal = 1, variability = [-10,10]%, 1 occurrences
  k: Uncertain real, nominal = 1, variability = [-10,10]%, 1 occurrences
  m: Uncertain real, nominal = 1, variability = [-10,10]%, 3 occurrences

Type "plant.NominalValue" to see the nominal value and "plant.Uncertainty" to interact with the uncertain elements.

From here, one can create the augmented plant (don't forget the frequency dependent weights!) for the D-K iterations.

[Peter DiFonzo]

Thank you so much!

I corrected the LFT matrix errors. I still get the same output, Krob = [ ].

I’ll try your suggestion. However, I still need to construct the basic mu synthesis architecture, (at the bottom of my write-up), correct? (I was following the procedure outlined in the textbook (and other sources)).

The application of frequency dependent weights is something I’ve been struggling with, however I gathered that I could add those later, once I understand the basic process.

[Paul]

Yes, you need to construct the basic mu-synthesis architecture as expected by musyn.

I just noticed another problem in the code ....

This line has the inputs backwards. The exogenous inputs in w should precede the control inputs in u in the call to connect, as shown in your block diagrams

P2=connect(a1,a2,Mm_sys,Mk_sys,Mc_sys,Sum1,Sum2,ap1,ap2,delta_sys,Sum3,{'u','r'},{'ye','u','ye','yx2'});

So it should be:

P2=connect(a1,a2,Mm_sys,Mk_sys,Mc_sys,Sum1,Sum2,ap1,ap2,delta_sys,Sum3,{'r','u'},{'ye','u','ye','yx2'});
%                                                                        ^   ^

The frequency dependent weights are going to be critical to get an acceptable design.

Feel free to post your updated code (don't need all of the graphics if they are too painful to add) in a comment response if you're still having trouble.

[Peter DiFonzo]

That did the trick!

Thanks so much!

I made the corrections in the code and the graphics, so others can follow along.

Now on to figuring out how to define and apply the weights, and figuring out how to interpret and implement the resultant Krob.

[Paul]

Hi Peter,

Based on the edits to the question it sounds like you've made significant progress.

I see that the edited code still doesn't address the first point in this answer, which is that the block diagram doesn't match the original differential equation.

From the question:

This system can be represented by the following differential equation of motion:

and represented by the following continuous time block diagram and state-space representation:

We see that the original differential equation is rewritten as

xddot = F/m - c/m*xdot - k/m*x

Therefore the k/m and c/m terms in the A-matrix should have negative signs and the summation in the block diagram to the right of the 1/m block should have negative feedback. If my assessment is correct, then you'd change

%Sum1=sumblk('s7=s9+u');
Sum1=sumblk('s7 = u - s9');

[Paul]

Unless there's a need to manually set up all of uncertainty blocks, I'd attack the problem this way (assuming standard mass-spring-damper model).

m0 = 1;
c0 = 1;
k0 = 1;
m = ureal("m",m0,'Percentage',10);
c = ureal("c",c0,'Percentage',10);
k = ureal("k",k0,'Percentage',10);
plant = ss([0 1;-k/m -c/m],[0;1/m],eye(2),0,'InputName','u','OutputName',{'x1','x2'});

errsum = sumblk('ye = r - x1');

Here is where the frequency dependent weights would be defined to defing the desired performance goals. Set to unity here to match your model up to this point

We = ss(1); We.InputName = 'ye';We.OutputName = 'yew';
Wu = ss(1); Wu.InputName = 'u'; Wu.OutputName = 'uw';

Connect the design plant. Include outputs for later analysis of the closed loop reponses to x1, ye, and u from r. The first three outputs are truth, the next two ouputs are the weighted signals in z, and the last two outputs are the measurements for feedback

sys = connect(plant,errsum,We,Wu,{'r','u'},{'x1', 'ye', 'u', 'yew', 'uw', 'ye','x2'});

Rename the measured error signal to avoid potential confusion.

sys.OutputName{end-1} = 'yemeas';

At this point, sys should be equivalent to P2 (if you change polarity on the feedback of Sum1) wrt to transfer functions from r and u to yew, uw, ye, and x2. Would be curious if that's actually the case.

Design the controller based only on the outputs in z

[K,CLperf,info] = musyn(sys({'yew', 'uw', 'yemeas','x2'},{'r','u'}),2,1);

D-K ITERATION SUMMARY:
-----------------------------------------------------------------
                       Robust performance               Fit order
-----------------------------------------------------------------
  Iter         K Step       Peak MU       D Fit             D
    1           1.003        1.001        1.001             0
    2           1.001            1            1             0
    3               1            1            1             0

Best achieved robust performance: 1

In principal, this controller should be identical to yours.

Form closed-loop system

syscl = lft(sys,K);

Step response

figure
step(syscl('x1','r'))

Frequency responses (only magnitude is important) from r to ye and u compared to inverses of the respective weights.

figure
bodemag(syscl('ye','r'),inv(We))

figure;
bodemag(syscl('u','r'),inv(Wu))

These plots would be much more interesting with appropriate definition of We and Wu.

[Peter DiFonzo]

Hi Paul,

Thank you so much the detailed explanation. Your methodology certainly looks like a more efficient way of approaching this problem. I'll work through your description.

Your help has been invaluable in helping me understand how to set up the musyn architecture! I’m eager to see if I get similar results.

-Pete

[Paul]

Hi Peter,

The current post states:

"Where Delta will be set to vary between -1 to +1,"

Based on that statement, and the valuse of the etas, I think the intent is to have the uncertain variables vary by +-10%. But these lines

delta_m = ureal("delta_m",1,'PlusMinus',1);% Uncertain LTI dynamics, +/-1
delta_k = ureal("deltal_k",1,'PlusMinus',1);% Uncertain LTI dynamics, +/-1 
delta_c = ureal("delta_c",1,'PlusMinus',1);% Uncertain LTI dynamics, +/-1 

set the range of each delta to 0 < delta < 2, which will make the uncertain parameters vary from 0 - 20% and the nominal parameters will 10% larger than they should be.

One potential advantage of your approach is that the uncertainty (delta_m) in the mass shows up as a single, scalar occurence in the model, whereas with my approach it's moded as a repeated scalar with multiplicity three. I think I've seen somewhere in the doc that the former may be preferable from a computational efficiency perspective.