home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
OS/2 Shareware BBS: 10 Tools
/
10-Tools.zip
/
octa21eb.zip
/
octave
/
SCRIPTS.ZIP
/
scripts
/
control
/
bddemo.m
< prev
next >
Wrap
Text File
|
1999-03-05
|
24KB
|
614 lines
# Copyright (C) 1996 A. Scottedward Hodel
#
# This file is part of Octave.
#
# Octave is free software; you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by the
# Free Software Foundation; either version 2, or (at your option) any
# later version.
#
# Octave is distributed in the hope that it will be useful, but WITHOUT
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
# for more details.
#
# You should have received a copy of the GNU General Public License
# along with Octave; see the file COPYING. If not, write to the Free
# Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
function bddemo()
# Octave Controls toolbox demo: Block Diagram Manipulations demo
# Written by David Clem August 15, 1994
# Modified by A S Hodel Summer-Fall 1996
str_sav = implicit_str_to_num_ok;
sav_page = page_screen_output;
implicit_str_to_num_ok = 1;
page_screen_output = 1;
while (1)
clc
k=0;
while(k > 14 || k < 1)
k = menu("Octave Block Diagram Manipulations Demo", ...
"sysadd/syssub: F(s) = G(s) +/- H(s)", ...
"sysappend: add new inputs/outputs", ...
"syssetsignals: change names of inputs, outputs, and/or states", ...
"sysconnect: connect specified system inputs/outputs", ...
"syscont/sysdisc: extract the continuous (discrete) part of a system", ...
"sysdup: duplicate specified inputs/outputs", ...
"sysgroup: group two systems into a single system,", ...
"sysmult: F(s) = G(s)*H(s) (series connection)", ...
"sysprune: keep only specified inputs/outputs", ...
"sysscale: scale inputs/outputs by specified gain matrices", ...
"parallel: parallel connection of two systems", ...
"buildssic: the combination of all", ...
"Design examples:", ...
"Return to main demo menu");
endwhile
if (k == 1)
clc
disp("sysadd: add two systems together")
disp("syssub: subtract F = G - H")
prompt
help sysadd
prompt
help syssub
prompt
disp('Example #1, \n')
cmd = "sys1 = tf2sys([1 -1],[1 2 1]);";
run_cmd
cmd = "sys2 = tf2sys([1 -1],[1 2 3]);";
run_cmd
disp("sys1=")
sysout(sys1);
prompt
disp("sys2=")
sysout(sys2);
cmd = "sys_sum1 = sysadd(sys1,sys1);";
run_cmd
disp("This example adds sys1 to itself.")
cmd = "sysout(sys_sum1)";
run_cmd
disp("Notice that the numerator is twice what it used to be.")
prompt
disp("Example 2:")
cmd = "sys_sub1 = syssub(sys1,sys2);";
run_cmd
disp("Notice that sysadd (actually sysgroup, called by sysadd) lets you")
disp("know that your two systems had identical names for their states,")
disp("inputs and outputs. The warning message is perfectly harmless,")
disp("and is provided for user information only.")
disp("sys_sub1=")
sysout(sys_sub1)
disp("Notice that, since the two transfer functions had different")
disp("denominators, syssub converted the primary system type to ")
disp("state space. This is the typical behavior unless the")
disp("the two systems are both SISO, they both have the same poles,")
disp("and at least one of them has tf for its primary system type.");
prompt
elseif (k == 2)
disp("sysappend: add new inputs and/or outputs to a system")
help sysappend
prompt
disp("Consider a double-integrator system:")
sys = tf2sys(1,[1 0 0]);
sys=sysupdate(sys,"ss");
sysout(sys,"ss");
disp("We add a velocity disturbance input as follows:")
cmd = "sys1=sysappend(sys,[1;0]);";
run_cmd
sysout(sys1,"ss");
disp("Names of inputs can be included as follows:")
cmd = "sys1=sysappend(sys,[1;0], [],[],[],\"Disturb\");";
run_cmd
disp("Notice that empty matrices can be listed for the D matrix if")
disp("all entries are zeros.")
disp(" ")
disp("sys1 is thus:")
sysout(sys1);
prompt
elseif (k == 3)
disp("syssetsignals:")
help syssetsignals
disp("Example system");
a = rand(3,3);
b = rand(3,2);
c = rand(2,3);
sys = ss2sys(a,b,c);
sysout(sys);
prompt
disp("Change state names to larry, moe, and curly as follows:")
cmd = "sys = syssetsignals(sys,\"st\",list(\"larry\",\"moe \" , \"curly\"));";
run_cmd
disp("Indicate that output 2 is discrete-time:")
cmd = "sys = syssetsignals(sys,\"yd\",1,2);";
run_cmd
disp("Change output 2 name to \"Vir\"");
cmd = "sys = syssetsignals(sys,\"out\",\"Vir\",2);";
run_cmd
disp("Resulting system is:")
sysout(sys);
prompt
elseif (k == 4)
help sysconnect
prompt
disp("********* N O T E *********")
disp("sysconnect is demonstrated fully in the design examples (option 13)");
prompt
elseif (k == 5)
disp("syscont and sysdisc: ")
disp("Example block diagram 1:")
disp(" ------------------ ---------------------");
disp(" u_in ->| Discrete system |--->| Continuous system | ---> y_out");
disp(" ------------------ ---------------------");
sys1 = tf2sys([1 2],[1 2 1], 1,"u_in","y_disc");
sys2 = tf2sys([1 0],[1 -3 -2],0,"c_in","y_out");
sys = sysmult(sys2,sys1);
disp("Consider the hybrid system")
sysout(sys);
prompt
help syscont
disp("The continuous part of the system can be extracted with syscont")
cmd = "[csys,Acd,Ccd] = syscont(sys);";
run_cmd
disp("The resulting csys is")
sysout(csys);
disp("Notice that B is 0; there is no purely continuous path from the")
disp("input to the output");
prompt
help sysdisc
disp("The discrete part of the system can be extracted with sysdisc")
cmd = "[dsys,Adc,Cdc] = sysdisc(sys)";
run_cmd
disp("The resulting dsys is")
sysout(dsys);
disp("sysdisc returns dsys=empty since sys has no discrete outputs.");
prompt
disp("Example block diagram 2:")
sys1 = tf2sys([1 2],[1 2 1], 1,"u_in","y_disc");
sys2 = tf2sys([1 0],[1 -3 -2],0,"c_in","y_out");
disp(" ---------------------")
disp(" u_in -->o-->| Discrete system | --------> y_disc")
disp(" ^ --------------------- |")
disp(" | | ");
disp(" -----------------------------|---")
disp(" | |")
disp(" ------------------------------ |")
disp(" | |")
disp(" v --------------------- |")
disp(" c_in -->o-->| continuous system | --------> y_out")
disp(" ---------------------")
disp("repeat the above example with sys=")
sys = sysgroup(sys1, sys2);
sysout(sys)
prompt
sys = sysconnect(sys,[1 2],[2 1]);
sysout(sys);
cmd = "[csys,Acd,Bcd] = syscont(sys);";
run_cmd
cmd = "[dsys,Acd,Bcd] = sysdisc(sys);";
run_cmd
disp("csys is now")
sysout(csys)
disp("dsys is now")
sysout(dsys);
prompt
elseif (k == 6)
help sysdup
prompt
disp("********* N O T E *********")
disp("sysdup is fully demonstrated in the design examples (option 13)")
prompt
elseif (k == 7)
help sysgroup
disp(" ")
prompt
disp("Example: combine two SISO systems together:")
cmd = "sys_a=tf2sys([1 2],[3 4]);";
run_cmd
cmd = "sys_b=tf2sys([5 6],[7 8],1);";
run_cmd
cmd = "sys_g=sysgroup(sys_a,sys_b);";
run_cmd
disp("Notice that sysgroup warns you when you join a purely continuous")
disp("system to a purely discrete system. sysgroup also warns when")
disp("you join two systems that have common state, input, or output names.")
cmd = "sysout(sys_g)";
run_cmd
disp("Since the grouped system is a multiple-input multiple-output system,")
disp("the output system is by default in state-space format.")
disp(" ")
disp("********* N O T E *********")
disp("sysgroup is further demonstrated in the design examples (option 13)")
prompt
elseif (k == 8)
help sysmult
disp("sysmult performs a series connection of two systems.")
disp("Example 1")
disp(" ")
disp(" ---------- ----------")
disp(" u --->| Bsys |---->| Asys |---> y")
disp(" ---------- ----------")
disp(" ")
Asys = tf2sys(1,[1 2 1],0,"a_in","a_out");
Bsys = tf2sys([2 3],[1 3 2],0,"b_in","b_out");
disp("Asys=")
sysout(Asys);
disp("Bsys=");
sysout(Bsys);
cmd = "sys = sysmult(Asys,Bsys);";
run_cmd
disp("sys =")
sysout(sys);
disp("Notice that sysmult automatically transforms to state space")
disp("internal representation. This is to avoid numerical problems")
disp("when multiplying polynomials");
prompt
disp("Example 2: same system, except that Bsys is discrete-time");
Bsys = tf2sys([2 3],[1 3 2],1e-2,"b_in","b_out");
sysout(Bsys);
cmd = "sys = sysmult(Asys,Bsys);";
run_cmd
disp("sys =")
sysout(sys);
prompt
elseif (k == 9)
help sysprune
prompt
disp("********* N O T E *********")
disp("sysprune is demonstrated in the design examples (option 13)");
prompt
elseif (k == 10)
help sysscale
disp(" ")
prompt
disp("********* N O T E *********")
disp("See the design examples (option 13) for use of sysscale.")
prompt
elseif ( k == 11)
help parallel
disp("parallel operates by making a call to sysgroup and sysscale.")
disp("Example:")
sys1 = tf2sys(1,[1 1],0,"in1","out1");
sys2 = tf2sys(2,[1 2],0,"in2","out2");
disp("sys1=")
sysout(sys1);
disp("sys2=")
sysout(sys2);
cmd = "sysp = parallel(sys1,sys2);";
run_cmd
disp("sysp=")
sysout(sysp);
prompt
disp("parallel can be used for multiple input systems as well:")
in1 = list("u1.1","u1.2");
in2 = list("u2.1","u2.2");
out1 = list("y1.1","y1.2");
out2 = list("y2.1","y2.2");
sys1 = ss2sys([-1,0;0 -2],eye(2),eye(2),[]);
sys2 = ss2sys([-2,0;0 -4],eye(2),eye(2),[]);
sys1 = syssetsignals(sys1,"in",in1);
sys1 = syssetsignals(sys1,"out",out1);
sys2 = syssetsignals(sys2,"in",in2);
sys2 = syssetsignals(sys2,"out",out2);
disp("sys1=")
sysout(sys1);
disp("sys2=")
sysout(sys2);
cmd = "sysp = parallel(sys1,sys2);";
run_cmd
disp("sysp=")
sysout(sysp);
prompt
elseif (k == 12)
# buildssic description
disp(" ")
disp(" ---------------------------------------")
disp(" b u i l d s s i c")
disp(" (BUILD State Space InterConnections)")
disp(" ---------------------------------------")
disp(" ")
disp("buildssic builds a single system from up to 8 systems.")
disp("It's primary pupose is the forming of interconnections")
disp("for H2/H_inf designs and the building of closed loop")
disp("systems.")
disp("The interconnections may be of arbitrary complexity.")
disp("The use of buildssic is an alternative to sysgroup,")
disp("sysadd/syssub, sysappend, sysconnect, sysdup, sysmult")
disp("sysprune, sysscale, parallel etc.")
disp("In contrast to these functions buildssic does not")
disp("handle mixed continuous and discrete systems. However,")
disp("discrete systems can be connected as long as their")
disp("sampling times are identical. Another drawback: the")
disp("names of input/output and state variables are clobbered.")
disp("Of course, buildsysic is useful in combination with sysgroup,")
disp("sysmult etc.")
prompt
disp("********* N O T E *********")
disp("buildssic is demonstrated in the design examples (option 13)");
prompt
elseif (k == 13)
disp("Design examples")
disp("Packed system matrices may be connected and manipulated")
disp("With the functions listed below:")
disp(" sysdup: duplicate selected inputs/outputs")
disp(" sysconnect: connect selected inputs/outputs")
disp(" sysgroup: group two systems together")
disp(" sysprune: prune a system to keep only selected inputs and outputs")
disp(" sysscale:pre/post multiply a system by constant matrices")
disp(" buildssic: connect systems with arbitrary complexity.")
prompt
disp("As a simple example, we will construct the system block ")
disp("diagram shown below ")
disp(" ")
disp(" + -------- --------");
disp(" r(t) ---> (+) --->| K(s) |--->| P(s) | ----> y(t)");
disp(" -^ -------- -------- |");
disp(" | |");
disp(" ------------------------------");
disp(" ")
disp("where P(s) is the plant, K(s) is the controller.")
prompt
disp("Simple example: P(s) is a first order lag, K(s) is a PI ")
disp("controller")
nump = 1;
denp = [1 1];
disp("P(s)=")
tfout(nump,denp)
numk = [1 1];
denk = [1 0];
disp("\nK(s)=")
tfout(numk,denk);
prompt
disp("We'll show three approaches. ")
P = tf2sys(nump,denp,0,"plant input","plant output");
K = tf2sys(numk, denk,0,"controller input","controller output");
meth = 0;
while(meth != 5)
disp("The first method consists of the following steps:")
disp(" step 1: create systems P and K")
disp(" step 2: group P and K together")
disp(" step 3: create a summing junction")
disp(" step 4: connect outputs to respective inputs")
disp(" step 5: prune the desired i/o connections")
disp("The second method is done as follows:")
disp(" step 1: create systems P and K and a summing block S")
disp(" step 2: connect P, K, and S in series")
disp(" step 3: connect y to inverted summing connection")
disp(" step 4: prune the desired i/o connections")
disp("The third method uses buildssic:")
disp(" step 1: GW = buildssic(...,K,P)")
disp(" ")
disp("Other design examples are in dgkfdemo (controldemo option 7)")
disp(" ")
meth = menu("Select design example method", ...
"Method 1 ", ...
"Method 1 (w/o algebraic loop warning)", ...
"Method 2", ...
"Method 3", ...
"Exit design examples");
if(meth == 1)
disp(" * * * Method 1 * * *")
disp(" ")
disp(" + -------- --------");
disp(" r(t) ---> (+) --->| K(s) |--->| P(s) | ----> y(t)");
disp(" -^ -------- -------- |");
disp(" | |");
disp(" ------------------------------");
disp(" ")
disp("Step 1: put plants in system format:");
nump
denp
cmd = "P = tf2sys(nump,denp,0,""plant input"",""plant output"");";
run_cmd
disp("P=")
sysout(P)
prompt
numk
denk
cmd = "K = tf2sys(numk, denk,0,""controller input"",""controller output"");";
run_cmd
sysout(K)
prompt
disp("Step 2: group the systems together")
cmd = "PK = sysgroup(P,K);";
run_cmd
disp("PK=")
sysout(PK);
prompt
disp(" ")
disp(" y2 u1")
disp(" + -------- --------");
disp(" r(t) ---> (+) --->| K(s) |--->| P(s) | ----> y(t)");
disp(" u2 -^ -------- -------- | y1");
disp(" u3 | |");
disp(" --------------------------------");
disp(" ")
disp("The controller has two inputs summed together, r(t)")
disp("and the negative feedback of y(t)")
disp("Step 3a: duplicate controller input: (input 2 of PK)")
prompt
cmd = "PK = sysdup(PK,[],2);";
run_cmd
disp("PK=")
sysout(PK);
disp("Notice that PK now has three inputs (input 3 is a duplicate ");
prompt("of input 2). Press return to go on")
disp("Step 3b: scale input 3 by -1")
cmd = "PK = sysscale(PK,[],diag([1,1,-1]));";
run_cmd
disp("PK=")
sysout(PK);
prompt
disp("Step 4: connect:")
disp(" y(t) (output 1) to the negative sum junction (input 3)")
disp(" u(t) (output 2) to plant input (input 1)")
disp("and prune extraneous inputs/outputs (retain input 2, output 1)")
prompt
out_connect = [1 2]
in_connect = [3 1]
cmd = "PK0 = sysconnect(PK,out_connect,in_connect);";
run_cmd
prompt
disp("Notice that sysconnect detects the possibility of algebraic")
disp("connections when connecting inputs. Option 2 (Method 1 ")
disp("without algebraic loops) shows how to avoid this warning")
disp("by performing connections one at a time.")
prompt
disp("PK0=")
sysout(PK0);
disp("Notice that the connected inputs now have stars on their")
disp("names. This is how the Octave controls toolbox reminds you")
disp("that the loop has been closed around these inputs.")
prompt("Press return to prune extra inputs/outputs from system")
disp("Only keep plant output (output 1) and r(t) (input 2)")
cmd = "PK0 = sysprune(PK0,1,2);";
run_cmd
disp("PK0=")
sysout(PK0);
prompt
disp("The resulting closed-loop transfer function is obtained as follows:")
cmd = "[num,den] = sys2tf(PK0);";
run_cmd
prompt
disp("Transfer function is now")
tfout(num,den)
disp("You can check this: Pk0=PK/(1+PK), as expected")
prompt
elseif(meth == 2)
disp("Method 1 without algebraic loops")
disp(" ")
disp(" y2 u1")
disp(" + -------- --------");
disp(" r(t) ---> (+) --->| K(s) |--->| P(s) | ----> y(t)");
disp(" u2 -^ -------- -------- | y1");
disp(" u3 | |");
disp(" --------------------------------");
disp(" ")
disp("Recall that sysconnect checks for algebraic loops. Although")
disp("Design option 1 gets a warning message about a possible");
disp("algebraic loop, such a loop does not exist.")
disp("This can be seen by performing the connections one at a time");
cmd = "PK = sysgroup(P,K);";
run_cmd
disp("PK=")
sysout(PK);
disp("Create an additial inverted input to the controller.")
cmd = "PK = sysdup(PK,[],2);";
run_cmd
cmd = "PK = sysscale(PK,[],diag([1,1,-1]));";
run_cmd
disp("PK=")
sysout(PK);
disp("Connect controller to plant:")
cmd = "PK0 = sysconnect(PK,2,1);";
run_cmd
disp("Plant output to negative control input")
cmd = "PK0 = sysconnect(PK0,1,3);";
run_cmd
disp("Only keep plant output (output 1) and r(t) (input 2)")
cmd = "PK0 = sysprune(PK0,1,2);";
run_cmd
disp("PK0=")
sysout(PK0);
prompt
disp("The transfer function form of PK0 is:")
sysout(PK0,"tf");
prompt
elseif(meth == 3)
disp(" * * * Method 2 * * *")
disp(" ")
disp(" + -------- --------");
disp(" r(t) ---> (+) --->| K(s) |--->| P(s) | ----> y(t)");
disp(" -^ -------- -------- |");
disp(" | |");
disp(" ------------------------------");
disp(" ")
disp("Step 1: We've already created systems P and K. Create a sum ")
disp("block as follows:")
implicit_str_to_num_ok = "warn";
cmd = "S = ss2sys([],[],[],[1 -1],0,0,0,[],list(""r(t)"",""y(t)""),""e(t)"");";
run_cmd
implicit_str_to_num_ok = 1;
disp("You may avoid the string conversion warning by setting the ")
disp("Octave global variables implicit_str_to_num_ok = 1");
disp(" ");
disp("(You may wish to look at help ss2sys to see what the above does)");
disp("S=")
sysout(S)
disp("notice that this is just a gain block that outputs e = r - y")
prompt
disp("Step 2: series connections of P, K, and S")
cmd = "PKS = sysmult(P,sysmult(K,S));";
run_cmd
disp("PKS=")
sysout(PKS)
disp("Step 3: connect y to inverted input")
cmd = "PKcl = sysconnect(PKS,1,2);";
run_cmd
disp("PKcl=")
sysout(PKcl)
disp("Step 4: prune desired inputs/outputs")
cmd = "PKcl=sysprune(PKcl,1,1);";
run_cmd
disp("PKcl=")
sysout(PKcl)
prompt
elseif(meth == 4)
disp(" * * * Method 3 * * *")
disp(" ")
disp(" + -------- --------");
disp(" r(t) ---> (+) --->| K(s) |--->| P(s) | ----> y(t)");
disp(" -^ -------- -------- |");
disp(" | |");
disp(" ------------------------------");
disp(" ")
disp("Step 1: We've already created systems P and K.")
disp(" Let us call buildssic:")
disp(" PKcl = buildssic([1 2;2 -1],[],[1],[2],P,K)")
disp(" ")
disp(" ^ ^ ^ ^ ^ ^")
disp(" | | | | | |")
disp(" Connection list ----+ | | | | |")
disp(" internal input list -----------+ | | | +-- controller")
disp(" output list --------------+ | |")
disp(" input list ------------------+ +---- plant")
disp(" ")
disp(" Connection list: connect input 1 (P) with output 2 (K)")
disp(" connect input 2 (K) with neg. outp. 1 (P)")
disp(" ")
disp(" int. inp. list: do not append internal inputs")
disp(" (e.g. the internal input of K (r-y))")
disp(" ")
disp(" output list: the only output is 1 (P), positive")
disp(" ")
disp(" input list: the only input is 2 (K), positive")
disp(" ")
cmd = "PKcl = buildssic([1 2;2 -1],[],[1],[2],P,K);"
run_cmd
sysout(PKcl)
prompt
disp("The transfer function form of PKcl is:")
sysout(PKcl,"tf");
disp("You can check this: PKcl = PK / (1 + PK), as expected")
prompt
elseif(meth != 5)
disp("Illegal selection")
endif
endwhile
elseif (k == 14)
return
endif
endwhile
implict_str_to_num_ok = str_sav;
page_screen_output = sav_page;
endfunction