using IPOPT: for the optimization. Its open-source and I dont need optimization toolbox anymore

parent 56ff22a4
......@@ -13,124 +13,143 @@
% der: derivate of the transmission with respcet to the coefficients in x
function [Filter,tx]=MaximizeTransmission(varargin)
% Default options
opt = optimset('GradObj','on','DerivativeCheck','on','MaxIter',5000, ...
'MaxFunEvals',2000,'Display','off','TolFun',1e-6,'TolX',1e-6);
p=inputParser;
p.addRequired( 'Object2', @(x) isa(x,'TwoPort'));
p.addRequired( 'InitialPoint', @(x) isa(x,'RationalFilter'));
p.addParameter('Stopband', [], @(x) isa(x,'FRM'));
p.addParameter('options', opt, @(x) isa(x,'struct'));
p.parse(varargin{:});
args = p.Results;
Filter = args.InitialPoint;
% args.Stopband = [];
% DissipationFactor
sigma = args.InitialPoint.Dissipation;
% Conpute passband (Frequency between -1 and 1)
ind = args.Object2.Freq_normalised<1 & args.Object2.Freq_normalised>-1;
Passband = args.Object2.Freq(ind);
% Get stopband
Stopband = args.Stopband.Freq;
% Get required rejection
Rejection = 10.^(double(args.Stopband.Data)./10);
Rejection = [zeros(1,length(Passband)),reshape(Rejection,1,[])];
% Resample Object 2 to evaluate in the stopband
args.Object2 = resample(args.Object2,[Passband(:); Stopband(:)]);
% Passband index
pass = [true(1,length(Passband)),false(1,length(Stopband))];
% Full frequency vector normalised (passband and stopband)
freq = reshape(args.Object2.Freq_normalised,1,[]);
% Target data
A11 = double(reshape(args.Object2.S.Data(1,1,:),1,[]));
A21 = double(reshape(args.Object2.S.Data(2,1,:),1,[]));
% Filter order
N=degree(args.InitialPoint.Q);
% Initial point
lambda = 0;
x0=[real(Filter.P), imag(Filter.P), abs(Filter.R(2:end)),lambda];
% Leading term of R
ldR = abs(Filter.R(1));
% Actual optimization
[x] = fminimax(@(x) fun(x),x0,[],[],[],[],[],[],...
@(x) transmission(x,ldR,N,A11,A21,freq,pass,sigma,Rejection ),args.options);
% Auxiliar data
auxdata.Freq = freq;
auxdata.A11 = A11;
auxdata.A21 = A21;
auxdata.N = N;
auxdata.LdR = ldR;
auxdata.Pass = pass;
auxdata.Sigma = sigma;
% Initial point: [ real(P), imag(P), R(2:end), lambda ]
x0=[real(Filter.P), imag(Filter.P), abs(Filter.R(2:end)),0];
% The callback functions.
funcs.objective = @objective;
funcs.gradient = @gradient;
funcs.constraints = @constraints;
funcs.jacobian = @jacobian;
funcs.jacobianstructure = @jacobianstructure;
% Bounds
options.lb = [-inf(1,length(x0)-1),0]; % Lower bound on the variables.
options.ub = [ inf(1,length(x0)-1),1]; % Upper bound on the variables.
options.cl = -inf(1,length(freq)); % Lower bounds on the constraints.
options.cu = Rejection; % Upper bounds on the constraints.
% Set the IPOPT options.
options.ipopt.hessian_approximation = 'limited-memory';
options.ipopt.mu_strategy = 'adaptive';
options.ipopt.tol = 1e-6;
options.auxdata = auxdata;
% Run IPOPT.
[x info] = ipopt_auxdata(x0,funcs,options)
% Transmission level
tx = x(end);
% Construct output object
p=x(1:N+1)+1i*x(N+2:2*N+2);
r=[double(Filter.R(1)), double(x(2*N+3:end-1).*exp(1i*angle(Filter.R(1))))];
tx = x(end);
Filter.P = Polynome(p);
Filter.R = Polynome(r);
Filter.Q = RationalFilter.polyq(Filter.P,Filter.R);
Filter.C=-Filter.R(1)./conj(Filter.R(1));
end
function [obj, der] = fun(x)
obj = -x(end);
der = [zeros(length(x)-1,1); -1];
% ----------------------------------------------------------------------
function f = objective (x,~)
f = -x(end);
end
function [Tx, eqc, derTx, eqder] = transmission(x,ldR,N,A11,A21,freq,pass,sigma,rejection)
% Degree of the transmission polynomial
nr=length(x)-2*N-3;
if nr<0
error('Vector of coefficients x must contain at least 2N+2 elements');
% ----------------------------------------------------------------------
function g = gradient (x,~)
g = [zeros(length(x)-1,1); -1];
end
% ----------------------------------------------------------------------
function c = constraints (x,aux)
% Construct reflection and transmission polynomials
p=Polynome(x(1:N+1)+1i*x(N+2:2*N+2));
r=Polynome([ldR, x(2*N+3:end-1)]);
p=Polynome(x(1:aux.N+1)+1i*x(aux.N+2:2*aux.N+2));
r=Polynome([aux.LdR, x(2*aux.N+3:end-1)]);
lambda = x(end);
Filter=RationalFilter(p,'R',r,'DissipationFactor',sigma);
Filter=RationalFilter(p,'R',r,'DissipationFactor',aux.Sigma);
% Filter response
S22=reshape(Filter.S22(freq),1,[]);
S21=reshape(Filter.S21(freq),1,[]);
S22=reshape(Filter.S22(aux.Freq),1,[]);
S21=reshape(Filter.S21(aux.Freq),1,[]);
% Compute transmission
Transmission=S21./(1-S22.*A11);
Tx1=(abs(Transmission).^2).*abs(A21).^2;
Transmission=S21./(1-S22.*aux.A11);
Tx=(abs(Transmission).^2).*abs(aux.A21).^2;
% Passband and stopband constraints
Tx = (lambda - Tx1).*pass + (Tx1 - rejection).*(~pass);
eqc = [];
% If nargout>2 compute the gradient of the transmission with respect of each
% coefficient (jacovian matrix)
if nargout>2
% Jacovian of parameters S22 and S21
J22=Filter.jacobianS22(freq);
J21=Filter.jacobianS21(freq);
% Derivative of the pseudohyperbolic distance
dT=(J21.*(1-S22.*A11)+J22.*S21.*A11)./((1-S22.*A11).^2);
derTx1=2.*real(conj(Transmission).*dT).*abs(A21).^2;
% Jacobian of passband and stopband constraints
derTx = -derTx1.*pass + Tx1.*(~pass);
% Derivative with respect to lambda
derTx = [derTx; ones(1,length(pass))];
eqder = [];
end
c = (lambda - Tx).*aux.Pass + Tx.*(~aux.Pass);
end
% ----------------------------------------------------------------------
function J = jacobian (x,aux)
% Construct reflection and transmission polynomials
p=Polynome(x(1:aux.N+1)+1i*x(aux.N+2:2*aux.N+2));
r=Polynome([aux.LdR, x(2*aux.N+3:end-1)]);
Filter=RationalFilter(p,'R',r,'DissipationFactor',aux.Sigma);
% Filter response
S22=reshape(Filter.S22(aux.Freq),1,[]);
S21=reshape(Filter.S21(aux.Freq),1,[]);
% Compute transmission
Transmission=S21./(1-S22.*aux.A11);
Tx=(abs(Transmission).^2).*abs(aux.A21).^2;
% Jacovian of parameters S22 and S21
J22=Filter.jacobianS22(aux.Freq);
J21=Filter.jacobianS21(aux.Freq);
% Derivative of the transmission
dTx=(J21.*(1-S22.*aux.A11)+J22.*S21.*aux.A11)./((1-S22.*aux.A11).^2);
dTx=2.*real(conj(Transmission).*dTx).*abs(aux.A21).^2;
% Jacobian of passband and stopband constraints
derTx = -dTx.*aux.Pass + dTx.*(~aux.Pass);
% Derivative with respect to lambda
derLambda = [ones(1,sum(aux.Pass)),zeros(1,sum(~aux.Pass))];
J = sparse([derTx; derLambda].');
end
function J = jacobianstructure(aux)
J = sparse(ones(length(aux.Freq),2*aux.N+3));
end
% Minimize the pseudohyperbolic distance between a filter and a TwoPort
%
% Input arguments:
% x: vector with the coefficients of the transmission (p) and reflection (r)
% polynomial of the filter in the format x=[real(p); imag(p); r(2:end)]
% ldR: leading term of the polynomial R
% N: degree of the polynomial p
% Target: Target function
% freq: frequency interval
%
% Output arguments:
% dist: pseudohyperbolic distance at the points in vector freq
% der: derivate of the distance with respcet to the coefficients in x
function [Filter,rx]=MinimizeReflection(varargin)
% Default options
opt = optimset('GradObj','on','DerivativeCheck','off','MaxIter',5000, ...
'MaxFunEvals',2000,'Display','off','TolFun',1e-6,'TolX',1e-6);
p=inputParser;
p.addRequired( 'Object2', @(x) isa(x,'TwoPort'));
p.addRequired( 'InitialPoint', @(x) isa(x,'RationalFilter'));
p.addParameter('options', opt, @(x) isa(x,'struct'));
p.parse(varargin{:});
args = p.Results;
Filter = args.InitialPoint;
% Get Frequency between -1 and 1
ind = args.Object2.Freq_normalised<1 & args.Object2.Freq_normalised>-1;
Freq = reshape(args.Object2.Freq_normalised(ind),1,[]);
freq = reshape(args.Object2.Freq_normalised(ind),1,[]);
% Target data
Target= reshape(conj(args.Object2.S.Data(1,1,ind)),1,[]);
target= reshape(conj(args.Object2.S.Data(1,1,ind)),1,[]);
% Filter order
N=degree(args.InitialPoint.Q);
% Leading term of R
ldR = abs(Filter.R(1));
% Auxiliar data
auxdata.Freq = freq;
auxdata.Target = target;
auxdata.N = N;
auxdata.LdR = ldR;
% Initial point: [ real(P), imag(P), R(2:end), lambda ]
x0=[real(Filter.P), imag(Filter.P), abs(Filter.R(2:end)),1];
% The callback functions.
funcs.objective = @objective;
funcs.gradient = @gradient;
funcs.constraints = @constraints;
funcs.jacobian = @jacobian;
funcs.jacobianstructure = @jacobianstructure;
% Bounds
options.lb = [-inf(1,length(x0)-1),0]; % Lower bound on the variables.
options.ub = [ inf(1,length(x0)-1),1]; % Upper bound on the variables.
options.cl = -inf(1,length(freq)); % Lower bounds on the constraints.
options.cu = zeros(1,length(freq)); % Upper bounds on the constraints.
% Set the IPOPT options.
options.ipopt.hessian_approximation = 'limited-memory';
options.ipopt.mu_strategy = 'adaptive';
options.ipopt.tol = 1e-6;
options.auxdata = auxdata;
% Run IPOPT.
[x info] = ipopt_auxdata(x0,funcs,options);
% Initial point
lambda = 1;
x0=[real(Filter.P), imag(Filter.P), abs(Filter.R(2:end)),lambda];
% Actual optimization
[x] = fminimax(@(x) fun(x),x0,[],[],[],[],[],[],...
@(x) pseudohyperbolic_distance(x,abs(Filter.R(1)),N,Target,Freq),args.options);
% Reflection level
rx = x(end);
% Construct output object
p=x(1:N+1)+1i*x(N+2:2*N+2);
r=[double(Filter.R(1)), double(x(2*N+3:end-1).*exp(1i*angle(Filter.R(1))))];
rx = x(end);
Filter.P = Polynome(p);
Filter.R = Polynome(r);
Filter.Q = RationalFilter.polyq(Filter.P,Filter.R);
Filter.C=-Filter.R(1)./conj(Filter.R(1));
end
function [obj, der] = fun(x)
obj = x(end);
der = [zeros(length(x)-1,1); 1];
% ----------------------------------------------------------------------
function f = objective (x,~)
f = x(end);
end
function [dist,eqc, der,eqder]=pseudohyperbolic_distance(x,ldR,N,Target,freq)
% Compute the pseudohyperbolic distance between a filter Target and a Target
%
% Input arguments:
% x: vector with the coefficients of the transmission (p) and reflection (r)
% polynomial of the filter in the format x=[real(p); imag(p); r(2:end)]
% ldR: leading term of the polynomial R
% N: degree of the polynomial p
% Target: Target function
% freq: frequency interval
%
% Output arguments:
% dist: pseudohyperbolic distance at the points in vector freq
% der: derivate of the distance with respcet to the coefficients in x
% Degree of the transmission polynomial
nr=length(x)-2*N-3;
if nr<0
error('Vector of coefficients x must contain at least 2N+2 elements');
% ----------------------------------------------------------------------
function g = gradient (x,~)
g = [zeros(length(x)-1,1); 1];
end
% ----------------------------------------------------------------------
function c = constraints (x,aux)
% Construct reflection and transmission polynomials
p=Polynome(x(1:N+1)+1i*x(N+2:2*N+2));
r=Polynome([ldR, x(2*N+3:end-1)]);
p=Polynome(x(1:aux.N+1)+1i*x(aux.N+2:2*aux.N+2));
r=Polynome([aux.LdR, x(2*aux.N+3:end-1)]);
lambda = x(end);
Filter=RationalFilter(p,'R',r);
% Filter response
S22=reshape(Filter.S22(freq),1,[]);
S22=reshape(Filter.S22(aux.Freq),1,[]);
% pseudohyperbolic distance
G=(S22-Target)./(1-S22.*conj(Target));
dist=abs(G).^2 - lambda;
eqc = [];
G=(S22-aux.Target)./(1-S22.*conj(aux.Target));
c=abs(G).^2 - lambda;
end
% ----------------------------------------------------------------------
function J = jacobian (x,aux)
% Construct reflection and transmission polynomials
p=Polynome(x(1:aux.N+1)+1i*x(aux.N+2:2*aux.N+2));
r=Polynome([aux.LdR, x(2*aux.N+3:end-1)]);
Filter=RationalFilter(p,'R',r);
% Filter response
S22=reshape(Filter.S22(aux.Freq),1,[]);
% If nargout=2 compute the gradient of the distance with respect of each
% coefficient (jacovian matrix)
% pseudohyperbolic distance
G=(S22-aux.Target)./(1-S22.*conj(aux.Target));
if nargout>2
% Jacovian of parameters S22
J=jacobianS22(Filter,freq);
% Jacovian of parameters S22
J=jacobianS22(Filter,aux.Freq);
% Derivative of the pseudohyperbolic distance
dG=J.*(1-abs(T).^2)./((1-S22.*conj(T)).^2);
der=2.*real(conj(G).*dG);
% Derivative of the pseudohyperbolic distance
dG=J.*(1-abs(aux.Target).^2)./((1-S22.*conj(aux.Target)).^2);
der=2.*real(conj(G).*dG);
% Derivative with respect to lambda
J = sparse([der; -ones(1,length(aux.Freq))].');
end
% Derivative with respect to lambda
der = [der; -ones(1,length(freq))];
eqder = [];
function J = jacobianstructure(aux)
J = sparse(ones(length(aux.Freq),2*aux.N+3));
end
end
\ No newline at end of file
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment