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

17 views (last 30 days)
Peter DiFonzo on 23 Apr 2024
Commented: Paul on 11 May 2024
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 on 11 May 2024
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.

Paul on 4 May 2024
Edited: Paul on 4 May 2024
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%.
Peter DiFonzo on 7 May 2024
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 on 7 May 2024
Edited: Paul on 7 May 2024
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.

### Categories

Find more on Uncertain Models in Help Center and File Exchange

R2023a

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!