function[fmin_str,fminnz_sol,clp_fpol]=sdsfminnz(E,A,b,c,clpf)
% This code will produce a gain vector fmin_str, which has minimun non-zero
% entries, for a given set of structured system matrices E,A,b.
% Further, for the exact numerical entries of E,A,b, and a set
% of desired closed loop finite poles, this code will also produce
% fminnz_sol, which is a numerical realization of fmin_str.
% The matrix c is required just to compute the finite poles of closed loop
% system using the command "dss".
% This version is modified on 19.12.2021 after 1st revision of Automatica
% review.
% This version is updated on 17.03.2021
% The code and comments of the code are updated on 21.09.2022. The new
% function "bptmaxmatch" is used to construct Zbt_str matrix, which
% utilizes the "maxflow" algorithm to compute a maximum matching of a
% bipartite graph.
%==========================================================================
sz_e = size(E);
E_sz = sz_e(1);
n = E_sz;
%% Extraction of Edge details of digraph G(H_0)
Eedg_s = [];
Eedg_t = [];
omega_e = [];
for i = 1:1:n
for j = 1:1:n
if (E(i,j) ~= 0) && (i ~= j)
Eedg_s = [Eedg_s j];
Eedg_t = [Eedg_t i];
omega_e = [omega_e E(i,j)];
end
end
end
lne = length(Eedg_s);
Eedg_ls = [];
Eedg_lt = [];
omega_el = [];
for i = 1:1:n
for j = 1:1:n
if (E(i,j) ~= 0) && (i == j)
Eedg_ls = [Eedg_ls j];
Eedg_lt = [Eedg_lt i];
omega_el = [omega_el E(i,j)];
end
end
end
lnel = length(Eedg_ls);
Aedg_s = [];
Aedg_t = [];
omega_a = [];
for i = 1:1:n
for j = 1:1:n
if (A(i,j) ~= 0) && (i ~= j)
Aedg_s = [Aedg_s j];
Aedg_t = [Aedg_t i];
omega_a = [omega_a -A(i,j)];
end
end
end
lna = length(Aedg_s);
Aedg_ls = [];
Aedg_lt = [];
omega_al = [];
for i = 1:1:n
for j = 1:1:n
if (A(i,j) ~= 0) && (i == j)
Aedg_ls = [Aedg_ls j];
Aedg_lt = [Aedg_lt i];
omega_al = [omega_al -A(i,j)];
end
end
end
lnal = length(Aedg_ls);
end_vtx = n+1;
Bedg_s = [];
Bedg_t = [];
omega_b = [];
for i = 1:1:n
if b(i,1) ~= 0
Bedg_s = [Bedg_s end_vtx];
Bedg_t = [Bedg_t i];
omega_b = [omega_b b(i,1)];
end
end
lnb = length(Bedg_s);
Fedg_s = 1:1:n;
Fedg_t = end_vtx*(ones(1,n));
omega_f = ones(1,n); % Set the edge weights for f-edge to 1.
lnf = length(Fedg_s);
e_s = [Eedg_s Aedg_s Bedg_s Fedg_s]; % Initial vertex indeces of the edges
e_t = [Eedg_t Aedg_t Bedg_t Fedg_t]; % Final vertex indeces of the edges
e_ls = [Eedg_ls Aedg_ls]; % Self-loop vertices
e_lt = [Eedg_lt Aedg_lt]; % Self-loop vertices
%% Construction of all spanning cycle families
[Wedg,Theta_sw,Theta_tw] = scfprledg(e_s,e_ls,e_t,e_lt)
Omega_edg = [omega_e omega_a omega_b omega_f omega_el omega_al]; % This
% gives the edge weights of all the edges in the digraph G(H_0).
Omega = [Wedg;Omega_edg];
str = sprank(E);
sz_theta = size(Theta_sw);
lthr = sz_theta(1);
lthc = sz_theta(2);
uid_e = Wedg(end,1:lne); % This gives the UINs assigned to E-edgs of G(H_0)
% Similarly, all edges UINs are extracted.
uid_a = Wedg(end,lne+1:lne+lna);
uid_b = Wedg(end,lne+lna+1:lne+lna+lnb);
uid_f = Wedg(end,lne+lna+lnb+1:lne+lna+lnb+lnf);
uid_el = Wedg(end,lne+lna+lnb+lnf+1:lne+lna+lnb+lnf+lnel);
uid_al = Wedg(end,lne+lna+lnb+lnf+lnel+1:lne+lna+lnb+lnf+lnel+lnal);
%% Construction of Zb and Z matrix from the SCFs, using Lemma 8
% Two matrices Zb and Z, having zero entries with size n*(r+1),
% are defined, and the entries are updated in the subsequent steps.
% To fill the entries of Zb, we need to know, which f-edge is used
% (let say f1i) and the number of E-edges are used (say k),
% then the nonzero entry for Zb is Zb(i,k+1). The matrix Z_01
Z = zeros(n,(str+1));
Zb = zeros(n,(str+1));
for i = 1:1:lthr
Thet_temp = [Theta_sw(i,:);Theta_tw(i,:)]; % ith SCF
Idthx = Theta_sw(i,end_vtx+1:end); % Obtained the unique
% idntification numbers of
% each edges in the ith SCF
num_e_edge = length(intersect([uid_e uid_el],Idthx)); % Gives the
% number of E-edges in ith SCF
int_wtf = intersect(Idthx,uid_f); % Gives the identification
% numbers of F-edges in ith SCF
zx_fp = find(Wedg(end,:) == int_wtf); % Gives the position of f-edges
% in Wedg (None of the SCFs can
% have more than one f-edge,
% since they all enter to the
% last vertex.)
f_indx = Wedg(1,zx_fp) ; % Gives the index of f-edge that
% is used in ith SCF
% The zero entries of Zb are now updated by adding 1 to its previous
% values, depending on the SCFs associated with it. For instance,
% Zb(2,3) = 2 implies there are two SCFs associated with the entry
% Zb(2,3). At the final step, we will get structured matrix Zb
% having zero-nonzero pattern.
Zb(f_indx,num_e_edge+1) = Zb(f_indx,num_e_edge+1)+1
% For exact numerical entries of the elements of the system structured
% matrices, a numerical realization matrix Z of structured matrix Zb
% is computed next.
[int_val,omg_po] = intersect(Omega(3,:),Idthx,'stable'); % Extracts the
% position of the edges in ith
% SCF in Omga matrix
omg_value = Omega(end,omg_po); % Extracts the weights of the
% edges associated with ith SCF
p_w = prod(omg_value); % product of the weights of the
% edges associated with ith SCF
scf_sv = Theta_sw(i,1:end_vtx);
scf_tv = Theta_tw(i,1:end_vtx);
num_vdcyl = numvdjcyl(scf_sv,scf_tv); % Number of vertex disjoint cycles
sign_wt = (-1)^(n+num_vdcyl) ;
p_sw = sign_wt*p_w;
Z(f_indx,num_e_edge+1) = Z(f_indx,num_e_edge+1)+p_sw % This
% gives the computed Z matrix
end
%% Computation of f^*, which has minimum free entries
syms 's'
al = det(s*E-A);
alpt = eval(coeffs(al,s));
str_deg_alp = length(alpt)-1;
alp = [];
Str_rankZb = sprank(Zb); % Gives structral rank of Z
if (Str_rankZb <= str)
disp(['Numerical realizations of SDS can NOT be made impulse-free and' ...
' place their finite poles at desired locations!'])
return
else
if (Str_rankZb == str+1) % for Not impulse-free and impulse-free systems
alp = [alpt zeros(1,str-str_deg_alp)];
[MM,Zb_str] = bptmaxmatch(Zb); % Structured Zb matrix which is used
% for f*min computation.
Zbt_str = Zb_str'
% % Zbt_str = strdrankmat(Zb) % Structured Z^T matrix which is used for
% f_min computation. This function is
% removed, and maxflow algorithm is used
% to compute maximum matching.
Zt = Z';
% Only the columns of Zt, which are associated with non-zero columns of
% Zbt_str, are considered and others are set to zeros in constructing
% Zt_str_fmin as follows. Similarly, fmin_str vector is also
% constructed. The 0-1 pattern of fmin_str is the
% structured gain vectior with minimum non-zero entries.
Zt_str_fmin = zeros((str+1),n);
fmin_str = zeros(1,n);
for j=1:1:E_sz
nz_Zbt_fmin = nnz(Zbt_str(:,j)) ;
if nz_Zbt_fmin ~= 0
Zt_str_fmin(:,j) = Zt(:,j)
fmin_str(1,j) = 1;
end
end
sig=fliplr(poly([clpf]));
as=(alp'-sig');
fminnz_sol = linsolve(Zt_str_fmin,as);
end
Gsc=dss((A+b*fminnz_sol'),b,c,0,E);
clp_fpol=pole(Gsc); % closed loop finite poles
end