function TxLine_Module1_Analytics_ML_v4_Polynomial
clc;
close all;
%% ============================================================
% TRANSMISSION LINE - MODULE 1
% MATLAB ONLINE COMPATIBLE VERSION
%
% Figures generated:
% 1. Voltage Waves
% 2. Current Waves
% 3. Standing Wave Magnitude
% 4. Power Along Line
% 5. ML Predicted vs Analytical
% 6. ML Accuracy By Parameter
%% ============================================================
%% ---------------- USER INPUTS ----------------
Z0_user = 50; % Characteristic impedance (ohm)
ZL_user = 100 + 1j*30; % Load impedance (ohm)
f_user = 1e9; % Frequency (Hz)
len_user = 0.10; % Physical line length (m)
vp_user = 2.0e8; % Phase velocity (m/s)
alpha_user = 0.20; % Attenuation constant (Np/m)
Zs_user = 50; % Source resistance (ohm)
Vs_user = 1; % Source RMS voltage (V)
%% ============================================================
% PART 1 - ANALYTICAL TRANSMISSION LINE CALCULATIONS
%% ============================================================
R = real(ZL_user);
X = imag(ZL_user);
beta_user = 2*pi*f_user/vp_user;
gamma_tl = alpha_user + 1j*beta_user;
% Reflection coefficient at load
GammaL = (ZL_user - Z0_user)/(ZL_user + Z0_user);
% Input impedance
tline = tanh(gamma_tl*len_user);
Zin = Z0_user * (ZL_user + Z0_user*tline) / ...
(Z0_user + ZL_user*tline);
% Output impedance
Zout = Z0_user * (Zs_user + Z0_user*tline) / ...
(Z0_user + Zs_user*tline);
% Reflection coefficient parameters
Gamma_mag = abs(GammaL);
Gamma_ang = angle(GammaL)*180/pi;
% VSWR
if Gamma_mag < 1
VSWR = (1 + Gamma_mag)/(1 - Gamma_mag);
else
VSWR = Inf;
end
% Return loss
if Gamma_mag > 0
ReturnLoss_dB = -20*log10(Gamma_mag);
else
ReturnLoss_dB = Inf;
end
%% ============================================================
% SOURCE / LOAD VOLTAGE AND CURRENT
%% ============================================================
Iin = Vs_user/(Zs_user + Zin);
Vin = Iin*Zin;
A = cosh(gamma_tl*len_user);
B = Z0_user*sinh(gamma_tl*len_user);
Vload = Vin/(A + B/ZL_user);
Iload = Vload/ZL_user;
Pin = real(Vin*conj(Iin));
Pload = real(Vload*conj(Iload));
if abs(Pin) > eps
efficiency = 100*Pload/Pin;
else
efficiency = NaN;
end
%% ============================================================
% DISPLAY NUMERICAL RESULTS
%% ============================================================
fprintf('\n');
fprintf('============================================================\n');
fprintf(' TRANSMISSION LINE MODULE-1 ANALYSIS\n');
fprintf('============================================================\n');
fprintf('Z0 = %.4f ohm\n',Z0_user);
fprintf('ZL = %.4f %+.4fj ohm\n', ...
real(ZL_user),imag(ZL_user));
fprintf('Frequency = %.4e Hz\n',f_user);
fprintf('Length = %.4f m\n',len_user);
fprintf('Phase velocity = %.4e m/s\n',vp_user);
fprintf('Alpha = %.4f Np/m\n',alpha_user);
fprintf('Beta = %.6f rad/m\n',beta_user);
fprintf('Gamma = %.6f %+.6fj\n', ...
real(GammaL),imag(GammaL));
fprintf('|Gamma| = %.6f\n',Gamma_mag);
fprintf('Gamma angle = %.4f deg\n',Gamma_ang);
fprintf('Zin = %.6f %+.6fj ohm\n', ...
real(Zin),imag(Zin));
fprintf('Zout = %.6f %+.6fj ohm\n', ...
real(Zout),imag(Zout));
fprintf('VSWR = %.6f\n',VSWR);
fprintf('Return Loss = %.6f dB\n',ReturnLoss_dB);
fprintf('Input Power = %.6e W\n',Pin);
fprintf('Load Power = %.6e W\n',Pload);
fprintf('Transmission Eff. = %.4f %%\n',efficiency);
fprintf('============================================================\n');
%% ============================================================
% PART 2 - WAVEFORM CALCULATIONS
%
% z = 0 at load
% z = length at source
%% ============================================================
z = linspace(0,len_user,1000);
Vplus_L = Vload/(1 + GammaL);
Vminus_L = GammaL*Vplus_L;
Vplus_z = Vplus_L .* exp(-gamma_tl*z);
Vminus_z = Vminus_L .* exp(+gamma_tl*z);
Vtotal_z = Vplus_z + Vminus_z;
Iplus_z = Vplus_z/Z0_user;
Iminus_z = -Vminus_z/Z0_user;
Itotal_z = Iplus_z + Iminus_z;
Pz = real(Vtotal_z .* conj(Itotal_z));
Vmag = abs(Vtotal_z);
%% ============================================================
% FIGURE 1 - VOLTAGE WAVES
%% ============================================================
figure('Name','Voltage Waves');
plot(z,abs(Vplus_z),'LineWidth',1.5);
hold on;
plot(z,abs(Vminus_z),'LineWidth',1.5);
plot(z,abs(Vtotal_z),'LineWidth',2);
grid on;
xlabel('Distance from Load (m)');
ylabel('|Voltage| (V)');
title('Transmission Line Voltage Waves');
legend('Incident Wave','Reflected Wave','Total Voltage', ...
'Location','best');
hold off;
%% ============================================================
% FIGURE 2 - CURRENT WAVES
%% ============================================================
figure('Name','Current Waves');
plot(z,abs(Iplus_z),'LineWidth',1.5);
hold on;
plot(z,abs(Iminus_z),'LineWidth',1.5);
plot(z,abs(Itotal_z),'LineWidth',2);
grid on;
xlabel('Distance from Load (m)');
ylabel('|Current| (A)');
title('Transmission Line Current Waves');
legend('Incident Current','Reflected Current','Total Current', ...
'Location','best');
hold off;
%% ============================================================
% FIGURE 3 - STANDING WAVE
%% ============================================================
figure('Name','Standing Wave Magnitude');
plot(z,Vmag,'LineWidth',2);
grid on;
xlabel('Distance from Load (m)');
ylabel('|V(z)| (V)');
title(sprintf('Standing Wave Magnitude - VSWR = %.3f',VSWR));
%% ============================================================
% FIGURE 4 - POWER ALONG TRANSMISSION LINE
%% ============================================================
figure('Name','Power Along Line');
plot(z,Pz,'LineWidth',2);
grid on;
xlabel('Distance from Load (m)');
ylabel('Power (W)');
title(sprintf('Power Along Transmission Line - Efficiency = %.2f%%', ...
efficiency));
%% ============================================================
% PART 3 - MACHINE LEARNING DATASET
%
% Reduced dataset size compared with original version so that
% MATLAB Online executes much faster.
%% ============================================================
fprintf('\nGenerating ML dataset...\n');
rng(42);
N = 4000;
%% Random transmission-line parameters
Z0 = 45 + 35*rand(N,1);
ZLr = 35 + 105*rand(N,1);
ZLi = -30 + 60*rand(N,1);
f = 0.80e9 + 0.40e9*rand(N,1);
L = 0.005 + 0.045*rand(N,1);
vp = 1.90e8 + 0.20e8*rand(N,1);
alpha = 0.01 + 0.20*rand(N,1);
Zs = 45 + 25*rand(N,1);
Vs = 0.80 + 0.80*rand(N,1);
%% Calculate physics
ZL = ZLr + 1j*ZLi;
beta = 2*pi*f./vp;
gg = alpha + 1j*beta;
t = tanh(gg.*L);
G = (ZL-Z0)./(ZL+Z0);
Zin_d = Z0.*(ZL + Z0.*t)./(Z0 + ZL.*t);
Zout_d = Z0.*(Zs + Z0.*t)./(Z0 + Zs.*t);
Gm = abs(G);
VSWR_d = (1+Gm)./max(1-Gm,1e-8);
RL_d = -20*log10(max(Gm,1e-12));
%% Voltage / power
A_d = cosh(gg.*L);
B_d = Z0.*sinh(gg.*L);
Iin_d = Vs./(Zs+Zin_d);
Vin_d = Iin_d.*Zin_d;
Vload_d = Vin_d./(A_d + B_d./ZL);
Iload_d = Vload_d./ZL;
Pin_d = real(Vin_d.*conj(Iin_d));
Pload_d = real(Vload_d.*conj(Iload_d));
eta_d = 100*Pload_d./max(abs(Pin_d),1e-12);
%% Output matrix
Y = [real(G), imag(G), ...
real(Zin_d), imag(Zin_d), ...
real(Zout_d), imag(Zout_d), ...
VSWR_d, RL_d, Pin_d, Pload_d, eta_d];
%% Input matrix
X = [Z0,ZLr,ZLi,f/1e9,L,vp/1e8,alpha,Zs,Vs];
%% ============================================================
% PART 4 - TRAIN / TEST SPLIT
%% ============================================================
idx = randperm(N);
nTrain = round(0.80*N);
tr = idx(1:nTrain);
te = idx(nTrain+1:end);
Xtr = X(tr,:);
Ytr = Y(tr,:);
Xte = X(te,:);
Yte = Y(te,:);
%% ============================================================
% NORMALIZE INPUTS
%% ============================================================
muX = mean(Xtr,1);
sigX = std(Xtr,[],1);
sigX(sigX < eps) = 1;
Xs_tr = (Xtr-muX)./sigX;
Xs_te = (Xte-muX)./sigX;
%% ============================================================
% POLYNOMIAL FEATURES
%
% Degree 3 is used here instead of degree 4 to make the
% MATLAB Online version substantially faster.
%% ============================================================
degree = 3;
exponents = makeExponentMatrix(size(Xs_tr,2),degree);
Phi_tr = makePolynomialFeatures(Xs_tr,exponents);
Phi_te = makePolynomialFeatures(Xs_te,exponents);
%% ============================================================
% STANDARDIZE POLYNOMIAL FEATURES
%% ============================================================
muPhi = mean(Phi_tr,1);
sigPhi = std(Phi_tr,[],1);
sigPhi(sigPhi < 1e-10) = 1;
Ps_tr = (Phi_tr-muPhi)./sigPhi;
Ps_te = (Phi_te-muPhi)./sigPhi;
%% ============================================================
% DIRECT ML TARGETS
%
% ReGamma
% ImGamma
% ReZin
% ImZin
% ReZout
% ImZout
% Pin
% Pload
% Efficiency
%% ============================================================
directCols = [1 2 3 4 5 6 9 10 11];
muYT = mean(Ytr(:,directCols),1);
sigYT = std(Ytr(:,directCols),[],1);
sigYT(sigYT < eps) = 1;
YT = (Ytr(:,directCols)-muYT)./sigYT;
%% ============================================================
% RIDGE REGRESSION
%% ============================================================
lambda = 1e-3;
nFeat = size(Ps_tr,2);
RidgeMatrix = lambda*eye(nFeat);
Beta = (Ps_tr'*Ps_tr + RidgeMatrix) \ ...
(Ps_tr'*YT);
%% ============================================================
% PREDICTION
%% ============================================================
predTeS = Ps_te*Beta;
predTeDirect = predTeS.*sigYT + muYT;
YtePred = zeros(size(Yte));
YtePred(:,directCols) = predTeDirect;
%% Calculate VSWR and Return Loss from predicted Gamma
GtePred = YtePred(:,1) + 1j*YtePred(:,2);
gmTe = abs(GtePred);
YtePred(:,7) = (1+gmTe)./max(1-gmTe,1e-8);
YtePred(:,8) = -20*log10(max(gmTe,1e-12));
%% ============================================================
% PART 5 - ACCURACY
%% ============================================================
nOut = size(Y,2);
r2 = zeros(1,nOut);
mae = zeros(1,nOut);
rmse = zeros(1,nOut);
for k = 1:nOut
residual = Yte(:,k)-YtePred(:,k);
ssres = sum(residual.^2);
sstot = sum((Yte(:,k)-mean(Yte(:,k))).^2);
r2(k) = 1-ssres/max(sstot,eps);
mae(k) = mean(abs(residual));
rmse(k) = sqrt(mean(residual.^2));
end
bestAcc = 100*mean(max(min(r2,1),0));
%% ============================================================
% DISPLAY ML RESULTS
%% ============================================================
outputNames = { ...
'ReGamma', ...
'ImGamma', ...
'ReZin', ...
'ImZin', ...
'ReZout', ...
'ImZout', ...
'VSWR', ...
'ReturnLoss', ...
'Pin', ...
'Pload', ...
'Efficiency'};
fprintf('\n');
fprintf('============================================================\n');
fprintf(' POLYNOMIAL ML VALIDATION\n');
fprintf('============================================================\n');
fprintf('Total data points = %d\n',N);
fprintf('Training samples = %d\n',numel(tr));
fprintf('Test samples = %d\n',numel(te));
fprintf('Polynomial degree = %d\n',degree);
fprintf('Polynomial features = %d\n',size(exponents,1));
fprintf('Ridge lambda = %.1e\n',lambda);
fprintf('Aggregate Test R2 = %.3f %%\n',bestAcc);
fprintf('============================================================\n');
fprintf('\nPer-output performance:\n');
for k = 1:nOut
fprintf('%-14s R2=%8.4f Accuracy=%8.3f%% MAE=%g RMSE=%g\n', ...
outputNames{k}, ...
r2(k), ...
100*max(min(r2(k),1),0), ...
mae(k), ...
rmse(k));
end
if bestAcc >= 97
fprintf('\n*** 97%% PROJECT ACCURACY TARGET ACHIEVED ***\n');
else
fprintf('\n*** Test accuracy is below 97%% in this run. ***\n');
end
%% ============================================================
% FIGURE 5 - ML PREDICTED VS ANALYTICAL
%% ============================================================
figure('Name','ML Predicted vs Analytical');
plotNames = { ...
'Re(Gamma)', ...
'Im(Gamma)', ...
'Re(Zin)', ...
'Im(Zin)'};
for q = 1:4
subplot(2,2,q);
plot(Yte(:,q),YtePred(:,q),'.');
hold on;
mn = min(Yte(:,q));
mx = max(Yte(:,q));
plot([mn mx],[mn mx],'--','LineWidth',1.5);
grid on;
xlabel(['Analytical ',plotNames{q}]);
ylabel(['ML ',plotNames{q}]);
title(sprintf('%s: R^2 = %.4f', ...
plotNames{q},r2(q)));
hold off;
end
%% ============================================================
% FIGURE 6 - ML ACCURACY
%% ============================================================
figure('Name','ML Accuracy By Parameter');
bar(1:nOut,100*max(min(r2,1),0));
hold on;
yline(97,'--','97% target','LineWidth',1.5);
grid on;
xticks(1:nOut);
xticklabels(outputNames);
xtickangle(45);
ylabel('R^2 Accuracy (%)');
xlabel('Output Parameter');
title(sprintf('ML Test Accuracy - Aggregate = %.2f%%',bestAcc));
hold off;
%% ============================================================
% FINISHED
%% ============================================================
fprintf('\n');
fprintf('============================================================\n');
fprintf(' ALL FIGURES GENERATED SUCCESSFULLY.\n');
fprintf('============================================================\n');
end
%% ============================================================
% HELPER FUNCTION 1
% Generate polynomial exponent matrix
%% ============================================================
function E = makeExponentMatrix(nVar,maxDegree)
E = [];
current = zeros(1,nVar);
E = generateExpRecursive(E,current,1,maxDegree);
function out = generateExpRecursive(out,vec,pos,remaining)
if pos == nVar
for a = 0:remaining
vec(pos) = a;
if sum(vec) > 0
out = [out;vec]; %#ok<AGROW>
end
end
else
for a = 0:remaining
vec(pos) = a;
out = generateExpRecursive( ...
out,vec,pos+1,remaining-a);
end
end
end
end
%% ============================================================
% HELPER FUNCTION 2
% Generate polynomial feature matrix
%% ============================================================
function P = makePolynomialFeatures(X,E)
n = size(X,1);
m = size(E,1);
P = ones(n,m);
for j = 1:m
v = ones(n,1);
for d = 1:size(X,2)
p = E(j,d);
if p ~= 0
v = v.*(X(:,d).^p);
end
end
P(:,j) = v;
end
end
To embed this project on your website, copy the following code and paste it into your website's HTML: