【无人机辅助网络定位】窃听环境下基于多臂赌博机方法的无人机辅助网络定位与功率优化研究(Matlab代码实现)

   💥💥💞💞欢迎来到本博客❤️❤️💥💥

🏆博主优势:🌞🌞🌞博客内容尽量做到思维缜密,逻辑清晰,为了方便读者。

⛳️座右铭:行百里者,半于九十。

📋📋📋本文内容如下:🎁🎁🎁

 ⛳️赠与读者

👨‍💻做科研,涉及到一个深在的思想系统,需要科研者逻辑缜密,踏实认真,但是不能只是努力,很多时候借力比努力更重要,然后还要有仰望星空的创新点和启发点。建议读者按目录次序逐一浏览,免得骤然跌入幽暗的迷宫找不到来时的路,它不足为你揭示全部问题的答案,但若能解答你胸中升起的一朵朵疑云,也未尝不会酿成晚霞斑斓的别一番景致,万一它给你带来了一场精神世界的苦雨,那就借机洗刷一下原来存放在那儿的“躺平”上的尘埃吧。

     或许,雨过云收,神驰的天地更清朗.......🔎🔎🔎

💥1 概述

窃听环境下基于多臂赌博机方法的无人机辅助网络定位与功率优化研究

摘要:本文针对窃听环境下的无人机辅助网络定位与功率优化问题展开研究。在存在窃听威胁的场景中,无人机作为移动节点,其定位精度与功率分配策略对网络通信的安全性和有效性至关重要。通过引入多臂赌博机方法,平衡无人机在定位探索与功率优化利用间的关系,提出一种自适应的定位与功率优化算法。仿真结果表明,该算法能够显著提升无人机网络在窃听环境下的定位精度与通信安全性,为复杂环境下的无线通信提供了一种有效的解决方案。

关键词:无人机辅助网络;窃听环境;多臂赌博机方法;定位优化;功率优化

一、引言

随着无线通信技术的快速发展,无人机辅助网络因其高机动性和灵活部署能力,在军事、民用等领域得到广泛应用。然而,在复杂的通信环境中,窃听行为成为威胁网络安全的重要因素。窃听者通过非法手段获取通信信息,不仅侵犯用户隐私,还可能导致关键信息泄露,对网络安全构成严重威胁。

在无人机辅助网络中,定位精度与功率分配策略直接影响通信质量与安全性。准确的定位能够确保无人机与地面用户之间的有效通信,而合理的功率分配则可以在满足通信需求的同时,降低被窃听的风险。多臂赌博机方法作为一种强化学习技术,能够有效平衡探索与利用的关系,为解决窃听环境下的定位与功率优化问题提供了新的思路。

二、相关研究现状

(一)无人机辅助网络定位研究

目前,无人机辅助网络定位方法主要包括相对视觉定位(RVL)和绝对视觉定位(AVL)。RVL方法通过帧到帧的匹配实现定位,如视觉里程计(VO)和同时定位与地图构建(SLAM)。然而,RVL方法存在累计漂移问题,长距离航行时定位误差会不断累积,无法满足高精度定位需求。AVL方法则根据无人机当前的观测在全局地图中进行匹配定位,现有的地图主要包括卫星地图和自行采集的地图。AVL方法的主要难点在于减少全局地图匹配时的搜索空间以及克服相机型号变化、四季变化导致的外观不一致问题。

(二)窃听环境下的通信安全研究

在窃听环境下,保障通信安全是无线通信系统的重要目标。现有研究主要从物理层安全角度出发,通过设计安全传输方案来提高通信的保密性。例如,无人机作为移动中继可以通过调整飞行轨迹靠近源节点收集信息,并在转发信息时远离窃听者位置,增加路径损耗以降低窃听信道质量;无人机搭载干扰器作为友好干扰,接近窃听节点广播人工噪声,干扰窃听者,使其信道质量变差。

(三)多臂赌博机方法的应用研究

多臂赌博机(MAB)是一种简化版的强化学习问题,主要用于解决在不确定环境下如何平衡探索与利用的决策问题。在MAB问题中,存在多个摇臂,每个摇臂对应一个未知的奖励分布。决策者需要在有限的尝试次数内,通过不断拉动摇臂,估计每个摇臂的期望奖励,并选择最优的摇臂以获得最大累计奖励。MAB方法在无线通信、推荐系统等领域得到广泛应用,为解决资源分配、决策优化等问题提供了有效的工具。

三、窃听环境下的无人机辅助网络定位与功率优化问题建模

(一)系统模型

考虑一个无人机辅助的无线通信网络,其中包含一个无人机、多个地面用户和一个潜在的窃听者。无人机作为移动节点,配备有通信设备和定位传感器,能够在三维空间中飞行。地面用户分布在特定区域内,通过与无人机进行通信来实现数据传输。窃听者试图窃听无人机与地面用户之间的通信信息。

假设无人机在任务时间T内,按照一定的轨迹飞行。将任务时间离散化为N个等距时隙,每个时隙长度为δt。在第n个时隙,无人机的位置表示为v[n] = [x[n], y[n], z[n]]T ∈ R³。地面用户的位置为u_k = [x_k, y_k]T ∈ R²,k ∈ [1, K],其中K为地面用户的数量。窃听者的位置为e = [x_e, y_e]T ∈ R²。

(二)定位模型

无人机通过接收地面用户发送的信号来实现定位。假设信号传播模型为视距(LoS)和非视距(NLoS)混合模型,其中LoS链路和NLoS链路的路径损耗分别由不同的公式表示。无人机在第n个时隙接收到来自第k个地面用户的信号强度为RSS_k[n],根据信号强度和已知的信号传播模型,可以估计无人机与地面用户之间的距离。

为了实现准确的定位,需要解决两个关键问题:一是如何减少在全局地图中匹配时的搜索空间;二是如何克服环境变化导致的外观不一致问题。本文采用基于特征点匹配的方法,结合三维地图信息,将无人机观测到的图像与地图中的特征进行匹配,从而实现定位。

(三)功率优化模型

在窃听环境下,无人机的功率分配策略需要同时考虑通信质量和安全性。一方面,要确保无人机与地面用户之间的通信质量,满足一定的信噪比(SNR)要求;另一方面,要降低被窃听者窃听的风险,即减小窃听者处的信干噪比(SINR)。

设无人机的发射功率为P,在第n个时隙,无人机与第k个地面用户之间的信道功率增益为g_uk[n],与窃听者之间的信道功率增益为g_ue[n]。地面用户处的噪声功率为σ²,窃听者处的噪声功率也为σ²。则地面用户处的信噪比为SNR_k[n] = P * g_uk[n] / σ²,窃听者处的信干噪比为SINR_e[n] = P * g_ue[n] / (σ² + I),其中I为其他干扰信号的功率。

功率优化的目标是在满足地面用户通信质量要求的前提下,最小化窃听者处的信干噪比,即:

其中,γ_{th}为地面用户通信所需的最低信噪比,P_{max}为无人机的最大发射功率。

(四)多臂赌博机问题建模

将无人机的定位探索和功率优化利用问题建模为多臂赌博机问题。每个定位策略和功率分配方案对应一个摇臂,拉动摇臂获得的奖励与定位精度和通信安全性相关。定义奖励函数为:

其中,Accuracy[n]为第n个时隙的定位精度,α和β为权重系数,用于平衡定位精度和通信安全性的重要性。

多臂赌博机问题的目标是在有限的尝试次数内,通过不断选择不同的定位策略和功率分配方案,最大化累计奖励,即:

四、基于多臂赌博机方法的定位与功率优化算法设计

(一)算法总体框架

本文提出的基于多臂赌博机方法的定位与功率优化算法主要包括三个模块:定位策略选择模块、功率分配方案选择模块和奖励评估模块。定位策略选择模块根据多臂赌博机方法选择最优的定位策略,功率分配方案选择模块选择最优的功率分配方案,奖励评估模块根据定位精度和通信安全性计算当前选择的奖励,并反馈给定位策略选择模块和功率分配方案选择模块,用于调整后续的选择策略。

(二)定位策略选择

采用ε-贪婪策略进行定位策略选择。在每个时隙,以概率ε随机选择一个定位策略进行探索,以概率1 - ε选择当前估计奖励最高的定位策略进行利用。定位策略包括基于特征点匹配的定位策略、基于模板匹配的定位策略等。通过不断探索和利用不同的定位策略,逐步估计每个定位策略的期望奖励,并选择最优的定位策略。

(三)功率分配方案选择

同样采用ε-贪婪策略进行功率分配方案选择。在每个时隙,以概率ε随机选择一个功率分配方案进行探索,以概率1 - ε选择当前估计奖励最高的功率分配方案进行利用。功率分配方案包括均匀功率分配、注水功率分配等。通过不断探索和利用不同的功率分配方案,逐步估计每个功率分配方案的期望奖励,并选择最优的功率分配方案。

(四)奖励评估与反馈

奖励评估模块根据定位精度和通信安全性计算当前选择的奖励。定位精度可以通过与地面真实位置的误差来衡量,通信安全性可以通过窃听者处的信干噪比来衡量。将计算得到的奖励反馈给定位策略选择模块和功率分配方案选择模块,用于更新每个定位策略和功率分配方案的期望奖励估计。

五、仿真实验与结果分析

(一)仿真环境设置

仿真实验在一个三维城市环境中进行,城市中包含多个建筑物。无人机在高度为H的平面上飞行,飞行轨迹为圆形轨迹。地面用户均匀分布在城市区域内,窃听者位于城市边缘。信号传播模型采用LoS和NLoS混合模型,参数根据实际环境进行设置。

(二)实验结果对比

将本文提出的基于多臂赌博机方法的定位与功率优化算法与传统的定位与功率优化算法进行对比。传统算法采用固定的定位策略和功率分配方案,不进行探索和利用的平衡。

实验结果表明,本文提出的算法在定位精度和通信安全性方面均优于传统算法。在定位精度方面,本文算法的平均定位误差比传统算法降低了约20%;在通信安全性方面,本文算法的平均窃听信干噪比比传统算法降低了约15%。

(三)算法性能分析

分析ε值对算法性能的影响。当ε值较大时,算法更倾向于探索,能够发现更好的定位策略和功率分配方案,但可能会导致短期内的奖励较低;当ε值较小时,算法更倾向于利用,能够快速收敛到当前最优解,但可能会陷入局部最优。通过实验发现,当ε = 0.1时,算法能够在定位精度和通信安全性之间取得较好的平衡。

六、结论与展望

(一)研究结论

本文针对窃听环境下的无人机辅助网络定位与功率优化问题,引入多臂赌博机方法,提出了一种自适应的定位与功率优化算法。通过平衡定位探索与功率优化利用间的关系,该算法能够显著提升无人机网络在窃听环境下的定位精度与通信安全性。仿真实验结果表明,本文算法优于传统算法,为复杂环境下的无线通信提供了一种有效的解决方案。

(二)研究展望

未来的研究可以进一步拓展以下几个方面:一是考虑更复杂的通信环境和窃听模型,如多窃听者、动态窃听等;二是结合其他强化学习算法,如深度强化学习,提高算法的性能和适应性;三是将算法应用到实际的无人机辅助网络中,进行实地测试和验证。

📚2 运行结果

主函数代码:

clc, clear, close all

%% Secrecy Improvement Block Scheme

%   Simulation for performance on SIB scheme evaluation.
%
%   A SIB is a Secrecy Improvement Block that is made up of nRLB RL blocks
%   in which nRL Bandits iterations are performed for a given dAB estimate,
%   and a Positioning block where the UAVs are made to assume new
%   positions.
%
%   Positions are based in 3 variables:
%       Rj:     Orbit radius of both UAVs around Alice
%       hj:     Height of both UAVs
%       theta:  Opening angle of UAVs behind Alice
%
%   Optimization method:    Non-static k-Armed Bandits UCB
%
%   Performance Variables
%       E_SIB:  Total energy consumed by the UAV at each SIB. Comprised of
%               Rx energy, positioning energy and ACK Tx energy.
%       S_SIB:  Total secrecy obtained. Basically a time integral of the
%               active WSC level over the SIB time.
%
%   Control Variables
%       sigAB:  Uncertainty in Bob's position.
%       nRLB:   Number of RL blocks in a SIB.
%   
%   Constants
%       nSIB:   Number of SIBs to be simulated. Could be adapted over
%               convergence criteria.
%       E_RX:   Energy needed by UAV to receive positioning information.
%               May be related to Rx power and number of bits for the
%               instruction.
%       E_ACK:  Same as E_RX but for the ACK.
%       P_Mov:  Power of UAV to move at a given speed v_J

%% Multi Varying nRLB
clc, clear, close all

% Environment parameters (Urban)

phi = 9.61;         % Environmental constant
omega = 0.16 ;      % Environmental constant
alpha = 0.3;        % Ground Path Loss exponent
alpha_AG = 0.3;     % Air-to-Ground Path Loss Exponent
ne_LOS = 1.0;       % Air-to-Ground LOS attenuation
ne_NLOS = 20;       % Air-to-Ground NLOS attenuation
KbVal = 5;          % K parameter for Bob's channels
KeVal = 5;          % K parameter for Eve's channels


% Channel parameters

m = 3;              % Number of parallel channels
sigma = 1/sqrt(2);  % Noise std dev of each component (real and imag) of every parallel channel

Rs = 1;


channelParam =  [   phi,...
                    omega,...
                    alpha,...
                    alpha_AG,...
                    ne_LOS,...
                    ne_NLOS,...
                    Rs,...
                    KbVal,...
                    KeVal];

% Positions ***************************************************************

    % Alice
nA = 1;             %   Number of Alices
A = [0,0,0];        %   Position of Alice (zero point)
gammaA = 200;        %   Alice Tx SNR

% -------------------------------------------------------------------------

    % Bob
nB = 1;             %   Number of Bobs
dAB_R = 20;        	%   Real distance between A and B
sigAB = dAB_R/10;          %   Unreliability of B's position
B = [dAB_R,0,0];      %   Position of Bob

% -------------------------------------------------------------------------

    % Eve
nR = 100;           %   Number of radial points
nTheta = 360;       %   Number of angular points
nE = nR*nTheta;     %   Number of Eves

rLow = 0.1;         %   Lowest radius of Eve
rHigh = 50;         %   Highest radius of Eve
thetaLow = 0;       %   Lowest angle of Eve
thetaHigh = 2*pi;   %   Highest radius of Eve

rangeR = linspace(rLow,rHigh,nR);                       % Points in Radial dimension
rangeTheta = linspace(thetaLow,thetaHigh,nTheta);       % Points in Angular dimension

thetat  =  repmat(rangeTheta,1,nR);
rt  = (repmat(rangeR',1,nTheta)).';

E = [rt(:).*cos(thetat(:)), rt(:).*sin(thetat(:)) , zeros(nR*nTheta,1)];             % Eves' position (rectangle coordinates)

% -------------------------------------------------------------------------

    % UAV
Nu = 2;                             %   Number of simultaneous UAVs
nUAV = 2;

nAng = 90;                          %   Angle discretization level (opening angle)  -> Number of Angle Actions
nH = 30;                            %   Height discretization level for both UAVs   -> Number of Height Actions
nRj = 30;                           %   Orbit radius discretization level           -> Number of radius Actions

gammaJ = 0.1*gammaA/nUAV;                 %   UAVs Jamming SNR

angleUAV    = linspace(0,2*pi/(nUAV-1),nAng);       %   Possible angle actions (opening angles)
hUAV        = linspace(100,200,nH);                    %   Possible height actions
rjUAV       = linspace(1,rHigh,nRj);                %   Possible orbit radius actions

% -------------------------------------------------------------------------

% *************************************************************************


%   UAV technical parameters [dummy values, research due]

E_RX = 0.1;                 % Energy it takes for UAV to receive positional information [dummy value].
E_ACK = 0.1;                % Energy it takes for UAV to send ACK to Alice [dummy value].
P_Mov = 80;                 % Power it takes UAV to maneuver [dummy value].
v_J = 2;                    % Speed at which UAV moves [dummy value].

%   SIB k-Armed Bandits
nSIB = 20;                  % Number of SIBs to be considered

nRLB_V = [5 10 20];
NnRLB = length(nRLB_V);

nRL = 1;                    % Number of RL steps per estimate
epsilon = 0;                % Epsilon for epsilon-greedy choice
initWSC = 0;                % Optimistic initial action values
c = 0.3;                    % Exploration parameter for UCB
alpha = 0.5;                  % Step size (0: uniform average)

%   Performance Variables

nMC = 1e2;                              %   Number of Monte Carlo loops for performance evaluation
E_SIB = zeros(NnRLB,nSIB);                  %   Energy used in a given SIB
S_SIB = zeros(NnRLB,nSIB);                  %   Cummulative secrecy over a given SIB

MC_WSC = zeros(NnRLB,nMC,nSIB);               %   Stored WSC per SIB for all MC loops
MC_Ang = zeros(NnRLB,nMC,nSIB);               %   Stored Angles per SIB for all MC loops
MC_H = zeros(NnRLB,nMC,nSIB);                 %   Stored Heights per SIB for all MC loops
MC_Rj = zeros(NnRLB,nMC,nSIB);                %   Stored Orbit Radius per SIB for all MC loops

dtRL = 0.1*3;                             % DUMMY TIME, to get consistent results, not depending on the processor or machine
                                        % Time per RL iteration (the smallest RL iteration)

boolFade = 1;               % Use fading

% Exhaustive Search
tic
fprintf('Exhaustive search!\n')
WSC_ES = exhaustiveSearchWSC(A, E, nUAV, angleUAV, hUAV, rjUAV, dAB_R, gammaA, gammaJ, channelParam, boolFade );
fprintf('\n\n Exhaustive search ended!')
toc

for iRLB = 1:NnRLB
    nRLB = nRLB_V(iRLB);
    nLoops = nSIB*(nRLB+1);     % Number of loops for action choosing (positioning every nRLB+1 loops)
    for iMC = 1:nMC
        % Initialization
        WSCEst_Angle    = initWSC*ones(1,nAng);         %   Action value estimation vector for angle actions
        WSCEst_H        = initWSC*ones(1,nH);           %   Action value estimation vector for height actions
        WSCEst_Rj       = initWSC*ones(1,nRj);          %   Action value estimation vector for orbit radius actions
        WSCN_Angle      = zeros(1,nAng);                %   Vector to store angle action ocurrences
        WSCN_H          = zeros(1,nH);                  %   Vector to store height action ocurrences
        WSCN_Rj         = zeros(1,nRj);                 %   Vector to store orbit radius action ocurrences

        WSC_UAV = 0;                                    %   WSC initial value
        
        WSC_Ang = angleUAV(end);               %   Angle random initialization
        WSC_H   = hUAV(end);                   %   Height random initialization
        WSC_Rj  = rjUAV(end);                  %   Orbit Radius random initialization

        tic
        for i=1:nLoops
            if mod(i,nRLB+1)==0
                % This is a positioning block
                % No RL is performed, and Performance variables are taken at
                % this block only
                t1 = toc;                                   % Time from previous SIB end to this Positioning block
                iSIB = i/(nRLB+1);                          % Number of SIB

                % Obtain new position indeces
                [indAng, indH, indRj] = getNewPos(WSCEst_Angle, WSCEst_H, WSCEst_Rj);   

                % UAV Movement
                WSC_Rj_Old = WSC_Rj;                                                    
                dH      = abs( WSC_H - hUAV(indH) );                                    % Height movement distance
                dRj     = abs( WSC_Rj - rjUAV(indRj) );                                 % Orbit radius movement distance
                dAng    = min(WSC_Rj,WSC_Rj_Old)*abs( WSC_Ang - angleUAV(indAng) )/2;	% Angle movement distance (through lesser orbit)
                dtJ     = (dAng + dH + dRj)/v_J;                                        % Flight time

                % Update UAV positioning values
                WSC_Ang = angleUAV(indAng);        
                WSC_H   = hUAV(indH);
                WSC_Rj  = rjUAV(indRj);

                % Update UAV position
                UAVs = setNewPos_N(nUAV, WSC_Ang, WSC_H, WSC_Rj);

                % Compute secrecy performance metric
                S_SIB(iRLB,iSIB) = S_SIB(iRLB,iSIB) +  dtRL*nRLB*3*WSC_UAV/nMC;

                % Compute energy metric
                E_SIB(iRLB,iSIB) = E_SIB(iRLB,iSIB) +  (E_RX + E_ACK + dtJ*P_Mov)/nMC;

                % Update UAV WSC value
                [WSC_UAV, DeltaComp] = computeWSC_NUAV(A, E, UAVs, dAB_R, gammaA, gammaJ, channelParam, boolFade );
                WSC_UAV = WSC_UAV/(nTheta*nR);

                % Store control values per MC loop
                MC_WSC(iRLB,iMC,iSIB)    = WSC_UAV;
                MC_Ang(iRLB,iMC,iSIB)    = WSC_Ang;               %   Stored Angles per SIB for all MC loops
                MC_H(iRLB,iMC,iSIB)      = WSC_H;                 %   Stored Heights per SIB for all MC loops
                MC_Rj(iRLB,iMC,iSIB)     = WSC_Rj;                %   Stored Orbit Radius per SIB for all MC loops

                fprintf('nRLB: %i \t\t MCLoop: %i of %i \t\t Positioning Block: %i of %i\t\tTime: %.3f\t\t WSC: %.3f\n',nRLB,iMC,nMC,iSIB,nSIB, t1,WSC_UAV);
                tic
                continue;
            end
            % If it is not a Positioning Block, it's an RL block
            % A single RL iteration is performed by each RL variable in
            % succession (angle, height, orbit radius) for a single dAB
            % estimate

            % dAB estimation and parameter computation
            dAB = normrnd(dAB_R,sigAB);                             %   CSI estimate

            %   RL iteration for angle (greedy on h and Rj)
            [WSCEst_Angle, WSCN_Angle]  = computeRLBlock_UCB_N(WSCEst_Angle, WSCEst_H, WSCEst_Rj, WSCN_Angle, angleUAV, hUAV, rjUAV, ...
                                                A, E, dAB, gammaA, gammaJ, c, channelParam, i, 1, alpha, nUAV, boolFade);

            %   RL iteration for H (greedy on angle and Rj)
            [WSCEst_H, WSCN_H]          = computeRLBlock_UCB_N(WSCEst_Angle, WSCEst_H, WSCEst_Rj, WSCN_H, angleUAV, hUAV, rjUAV, ...
                                                A, E, dAB, gammaA, gammaJ, c, channelParam, i, 2, alpha, nUAV, boolFade);

            %   RL iteration for Rj (greedy on angle and H)
            [WSCEst_Rj, WSCN_Rj]        = computeRLBlock_UCB_N(WSCEst_Angle, WSCEst_H, WSCEst_Rj, WSCN_Rj, angleUAV, hUAV, rjUAV, ...
                                                A, E, dAB, gammaA, gammaJ, c, channelParam, i, 3, alpha, nUAV, boolFade);
        end
    end
    save(['SIB_EX_nRLB_GenConst-v20nNF_', num2str(iRLB), '.mat'])
end

save('SIB_EX_nRLB_nUAV2-v20-nofade.mat')


%% Multi Varying sigma
clc, clear, close all

% Environment parameters (Urban)

phi = 9.61;         % Environmental constant
omega = 0.16 ;      % Environmental constant
alpha = 0.3;        % Ground Path Loss exponent
alpha_AG = 0.3;     % Air-to-Ground Path Loss Exponent
ne_LOS = 1.0;       % Air-to-Ground LOS attenuation
ne_NLOS = 20;       % Air-to-Ground NLOS attenuation
KbVal = 5;          % K parameter for Bob's channels
KeVal = 5;          % K parameter for Eve's channels


% Channel parameters

m = 3;              % Number of parallel channels
sigma = 1/sqrt(2);  % Noise std dev of each component (real and imag) of every parallel channel

Rs = 1;


channelParam =  [   phi,...
                    omega,...
                    alpha,...
                    alpha_AG,...
                    ne_LOS,...
                    ne_NLOS,...
                    Rs,...
                    KbVal,...
                    KeVal];

% Positions ***************************************************************

    % Alice
nA = 1;             %   Number of Alices
A = [0,0,0];        %   Position of Alice (zero point)
gammaA = 200;        %   Alice Tx SNR

% -------------------------------------------------------------------------

    % Bob
nB = 1;             %   Number of Bobs
dAB_R = 20;        	%   Real distance between A and B

sigAB_V = [dAB_R/40 dAB_R/10 dAB_R/5];
nSig = length(sigAB_V);

B = [dAB_R,0,0];      %   Position of Bob

% -------------------------------------------------------------------------

    % Eve
nR = 100;           %   Number of radial points
nTheta = 360;       %   Number of angular points
nE = nR*nTheta;     %   Number of Eves

rLow = 0.1;         %   Lowest radius of Eve
rHigh = 50;         %   Highest radius of Eve
thetaLow = 0;       %   Lowest angle of Eve
thetaHigh = 2*pi;   %   Highest radius of Eve

rangeR = linspace(rLow,rHigh,nR);                       % Points in Radial dimension
rangeTheta = linspace(thetaLow,thetaHigh,nTheta);       % Points in Angular dimension

thetat  =  repmat(rangeTheta,1,nR);
rt  = (repmat(rangeR',1,nTheta)).';

E = [rt(:).*cos(thetat(:)), rt(:).*sin(thetat(:)) , zeros(nR*nTheta,1)];             % Eves' position (rectangle coordinates)

% -------------------------------------------------------------------------

    % UAV
Nu = 2;                             %   Number of simultaneous UAVs
nUAV = 2;

nAng = 90;                          %   Angle discretization level (opening angle)  -> Number of Angle Actions
nH = 30;                            %   Height discretization level for both UAVs   -> Number of Height Actions
nRj = 30;                           %   Orbit radius discretization level           -> Number of radius Actions


gammaJ = 0.1*gammaA/nUAV;                 %   UAVs Jamming SNR

angleUAV    = linspace(0,2*pi/(nUAV-1),nAng);       %   Possible angle actions (opening angles)
hUAV        = linspace(100,200,nH);                    %   Possible height actions
rjUAV       = linspace(1,rHigh,nRj);                %   Possible orbit radius actions

% -------------------------------------------------------------------------

% *************************************************************************


%   UAV technical parameters [dummy values, research due]

E_RX = 0.1;                 % Energy it takes for UAV to receive positional information [dummy value].
E_ACK = 0.1;                % Energy it takes for UAV to send ACK to Alice [dummy value].
P_Mov = 80;                 % Power it takes UAV to maneuver [dummy value].
v_J = 2;                    % Speed at which UAV moves [dummy value].

%   SIB k-Armed Bandits
nSIB = 20;                  % Number of SIBs to be considered
nRLB = 10; 
nLoops = nSIB*(nRLB+1);     % Number of loops for action choosing (positioning every nRLB+1 loops)

nRL = 1;                    % Number of RL steps per estimate
epsilon = 0;                % Epsilon for epsilon-greedy choice
initWSC = 0;                % Optimistic initial action values
c = 0.3;                    % Exploration parameter for UCB
alpha = 0.5;                  % Step size (0: uniform average)

%   Performance Variables

nMC = 1e2;                              %   Number of Monte Carlo loops for performance evaluation
E_SIB = zeros(nSig,nSIB);                  %   Energy used in a given SIB
S_SIB = zeros(nSig,nSIB);                  %   Cummulative secrecy over a given SIB

MC_WSC = zeros(nSig,nMC,nSIB);               %   Stored WSC per SIB for all MC loops
MC_Ang = zeros(nSig,nMC,nSIB);               %   Stored Angles per SIB for all MC loops
MC_H = zeros(nSig,nMC,nSIB);                 %   Stored Heights per SIB for all MC loops
MC_Rj = zeros(nSig,nMC,nSIB);                %   Stored Orbit Radius per SIB for all MC loops

dtRL = 0.1*3;                             % DUMMY TIME, to get consistent results, not depending on the processor or machine
                                        % Time per RL iteration (the smallest RL iteration)

boolFade = 1;               % Use fading

% Exhaustive Search
tic
fprintf('Exhaustive search!\n')
WSC_ES = exhaustiveSearchWSC(A, E, nUAV, angleUAV, hUAV, rjUAV, dAB_R, gammaA, gammaJ, channelParam, boolFade );
fprintf('\n\n Exhaustive search ended!')
toc

for iRLB = 1:nSig
    sigAB = sigAB_V(iRLB); 
    for iMC = 1:nMC
        % Initialization
        WSCEst_Angle    = initWSC*ones(1,nAng);         %   Action value estimation vector for angle actions
        WSCEst_H        = initWSC*ones(1,nH);           %   Action value estimation vector for height actions
        WSCEst_Rj       = initWSC*ones(1,nRj);          %   Action value estimation vector for orbit radius actions
        WSCN_Angle      = zeros(1,nAng);                %   Vector to store angle action ocurrences
        WSCN_H          = zeros(1,nH);                  %   Vector to store height action ocurrences
        WSCN_Rj         = zeros(1,nRj);                 %   Vector to store orbit radius action ocurrences

        WSC_UAV = 0;                                    %   WSC initial value
        
        WSC_Ang = angleUAV(end);               %   Angle random initialization
        WSC_H   = hUAV(end);                   %   Height random initialization
        WSC_Rj  = rjUAV(end);                  %   Orbit Radius random initialization

        tic
        for i=1:nLoops
            if mod(i,nRLB+1)==0
                % This is a positioning block
                % No RL is performed, and Performance variables are taken at
                % this block only
                t1 = toc;                                   % Time from previous SIB end to this Positioning block
                iSIB = i/(nRLB+1);                          % Number of SIB

                % Obtain new position indeces
                [indAng, indH, indRj] = getNewPos(WSCEst_Angle, WSCEst_H, WSCEst_Rj);   

                % UAV Movement
                WSC_Rj_Old = WSC_Rj;                                                    
                dH      = abs( WSC_H - hUAV(indH) );                                    % Height movement distance
                dRj     = abs( WSC_Rj - rjUAV(indRj) );                                 % Orbit radius movement distance
                dAng    = min(WSC_Rj,WSC_Rj_Old)*abs( WSC_Ang - angleUAV(indAng) )/2;	% Angle movement distance (through lesser orbit)
                dtJ     = (dAng + dH + dRj)/v_J;                                        % Flight time

                % Update UAV positioning values
                WSC_Ang = angleUAV(indAng);        
                WSC_H   = hUAV(indH);
                WSC_Rj  = rjUAV(indRj);

                % Update UAV position
                UAVs = setNewPos_N(nUAV, WSC_Ang, WSC_H, WSC_Rj);

                % Compute secrecy performance metric
                S_SIB(iRLB,iSIB) = S_SIB(iRLB,iSIB) +  dtRL*nRLB*3*WSC_UAV/nMC;

                % Compute energy metric
                E_SIB(iRLB,iSIB) = E_SIB(iRLB,iSIB) +  (E_RX + E_ACK + dtJ*P_Mov)/nMC;

                % Update UAV WSC value
                [WSC_UAV, DeltaComp] = computeWSC_NUAV(A, E, UAVs, dAB_R, gammaA, gammaJ, channelParam, boolFade );
                WSC_UAV = WSC_UAV/(nTheta*nR);

                % Store control values per MC loop
                MC_WSC(iRLB,iMC,iSIB)    = WSC_UAV;
                MC_Ang(iRLB,iMC,iSIB)    = WSC_Ang;               %   Stored Angles per SIB for all MC loops
                MC_H(iRLB,iMC,iSIB)      = WSC_H;                 %   Stored Heights per SIB for all MC loops
                MC_Rj(iRLB,iMC,iSIB)     = WSC_Rj;                %   Stored Orbit Radius per SIB for all MC loops

                fprintf('sigma: %.3f \t\t MCLoop: %i of %i \t\t Positioning Block: %i of %i\t\tTime: %.3f\t\t WSC: %.3f\n',sigAB,iMC,nMC,iSIB,nSIB, t1,WSC_UAV);
                tic
                continue;
            end
            % If it is not a Positioning Block, it's an RL block
            % A single RL iteration is performed by each RL variable in
            % succession (angle, height, orbit radius) for a single dAB
            % estimate

            % dAB estimation and parameter computation
            dAB = normrnd(dAB_R,sigAB);                             %   CSI estimate

            %   RL iteration for angle (greedy on h and Rj)
            [WSCEst_Angle, WSCN_Angle]  = computeRLBlock_UCB_N(WSCEst_Angle, WSCEst_H, WSCEst_Rj, WSCN_Angle, angleUAV, hUAV, rjUAV, ...
                                                A, E, dAB, gammaA, gammaJ, c, channelParam, i, 1, alpha, nUAV, boolFade);

            %   RL iteration for H (greedy on angle and Rj)
            [WSCEst_H, WSCN_H]          = computeRLBlock_UCB_N(WSCEst_Angle, WSCEst_H, WSCEst_Rj, WSCN_H, angleUAV, hUAV, rjUAV, ...
                                                A, E, dAB, gammaA, gammaJ, c, channelParam, i, 2, alpha, nUAV, boolFade);

            %   RL iteration for Rj (greedy on angle and H)
            [WSCEst_Rj, WSCN_Rj]        = computeRLBlock_UCB_N(WSCEst_Angle, WSCEst_H, WSCEst_Rj, WSCN_Rj, angleUAV, hUAV, rjUAV, ...
                                                A, E, dAB, gammaA, gammaJ, c, channelParam, i, 3, alpha, nUAV, boolFade);
        end
    end
    save(['SIB_EX_nSig_GenConst-v0F_', num2str(iRLB), '.mat'])
end

save('SIB_EX_nSig_nUAV2-v10-fade.mat')

%% Multi Varying K (same for B and E)
clc, clear, close all

% Environment parameters (Urban)

phi = 9.61;         % Environmental constant
omega = 0.16 ;      % Environmental constant
alpha = 0.3;        % Ground Path Loss exponent
alpha_AG = 0.3;     % Air-to-Ground Path Loss Exponent
ne_LOS = 1.0;       % Air-to-Ground LOS attenuation
ne_NLOS = 20;       % Air-to-Ground NLOS attenuation

KbVal_v = [0, 5, 10];
KeVal_v = [0, 5, 10];

nKval = length(KbVal_v);
% Channel parameters

m = 3;              % Number of parallel channels
sigma = 1/sqrt(2);  % Noise std dev of each component (real and imag) of every parallel channel

Rs = 1;


% Positions ***************************************************************

    % Alice
nA = 1;             %   Number of Alices
A = [0,0,0];        %   Position of Alice (zero point)
gammaA = 200;        %   Alice Tx SNR

% -------------------------------------------------------------------------

    % Bob
nB = 1;             %   Number of Bobs
dAB_R = 20;        	%   Real distance between A and B
sigAB = dAB_R/10;
B = [dAB_R,0,0];      %   Position of Bob

% -------------------------------------------------------------------------

    % Eve
nR = 100;           %   Number of radial points
nTheta = 360;       %   Number of angular points
nE = nR*nTheta;     %   Number of Eves

rLow = 0.1;         %   Lowest radius of Eve
rHigh = 50;         %   Highest radius of Eve
thetaLow = 0;       %   Lowest angle of Eve
thetaHigh = 2*pi;   %   Highest radius of Eve

rangeR = linspace(rLow,rHigh,nR);                       % Points in Radial dimension
rangeTheta = linspace(thetaLow,thetaHigh,nTheta);       % Points in Angular dimension

thetat  =  repmat(rangeTheta,1,nR);
rt  = (repmat(rangeR',1,nTheta)).';

E = [rt(:).*cos(thetat(:)), rt(:).*sin(thetat(:)) , zeros(nR*nTheta,1)];             % Eves' position (rectangle coordinates)

% -------------------------------------------------------------------------

    % UAV
Nu = 2;                             %   Number of simultaneous UAVs
nUAV = 2;

nAng = 90;                          %   Angle discretization level (opening angle)  -> Number of Angle Actions
nH = 30;                            %   Height discretization level for both UAVs   -> Number of Height Actions
nRj = 30;                           %   Orbit radius discretization level           -> Number of radius Actions

gammaJ = 0.1*gammaA/nUAV;                 %   UAVs Jamming SNR

angleUAV    = linspace(0,2*pi/(nUAV-1),nAng);       %   Possible angle actions (opening angles)
hUAV        = linspace(100,200,nH);                    %   Possible height actions
rjUAV       = linspace(1,rHigh,nRj);                %   Possible orbit radius actions

% -------------------------------------------------------------------------

% *************************************************************************


%   UAV technical parameters [dummy values, research due]

E_RX = 0.1;                 % Energy it takes for UAV to receive positional information [dummy value].
E_ACK = 0.1;                % Energy it takes for UAV to send ACK to Alice [dummy value].
P_Mov = 80;                 % Power it takes UAV to maneuver [dummy value].
v_J = 2;                    % Speed at which UAV moves [dummy value].

%   SIB k-Armed Bandits
nSIB = 20;                  % Number of SIBs to be considered
nRLB = 10; % ------------------------------------------------------------------------------------- CHECK THIS
nLoops = nSIB*(nRLB+1);     % Number of loops for action choosing (positioning every nRLB+1 loops)
nRL = 1;                    % Number of RL steps per estimate
epsilon = 0;                % Epsilon for epsilon-greedy choice
initWSC = 0;                % Optimistic initial action values
c = 0.3;                    % Exploration parameter for UCB
alpha = 0.5;                  % Step size (0: uniform average)

%   Performance Variables

nMC = 1e2;                              %   Number of Monte Carlo loops for performance evaluation
E_SIB = zeros(nKval,nSIB);                  %   Energy used in a given SIB
S_SIB = zeros(nKval,nSIB);                  %   Cummulative secrecy over a given SIB

MC_WSC = zeros(nKval,nMC,nSIB);               %   Stored WSC per SIB for all MC loops
MC_Ang = zeros(nKval,nMC,nSIB);               %   Stored Angles per SIB for all MC loops
MC_H = zeros(nKval,nMC,nSIB);                 %   Stored Heights per SIB for all MC loops
MC_Rj = zeros(nKval,nMC,nSIB);                %   Stored Orbit Radius per SIB for all MC loops

dtRL = 0.1*3;                             % DUMMY TIME, to get consistent results, not depending on the processor or machine
                                        % Time per RL iteration (the smallest RL iteration)

boolFade = 1;               % Use fading

WSC_ES = zeros(nKval, length(angleUAV), length(hUAV), length(rjUAV));

for iRLB = 1:nKval
    KbVal = KbVal_v(iRLB);          % K parameter for Bob's channels
    KeVal = KeVal_v(iRLB);          % K parameter for Eve's channels
    
    channelParam =  [   phi,...
                    omega,...
                    0.3,...
                    alpha_AG,...
                    ne_LOS,...
                    ne_NLOS,...
                    Rs,...
                    KbVal,...
                    KeVal];

    % Exhaustive Search
    tic
    fprintf('Exhaustive search!\n')
    WSC_ES(iRLB,:,:,:) = exhaustiveSearchWSC(A, E, nUAV, angleUAV, hUAV, rjUAV, dAB_R, gammaA, gammaJ, channelParam, boolFade );
    fprintf('\n\n Exhaustive search ended!')
    toc
    for iMC = 1:nMC
        % Initialization
        WSCEst_Angle    = initWSC*ones(1,nAng);         %   Action value estimation vector for angle actions
        WSCEst_H        = initWSC*ones(1,nH);           %   Action value estimation vector for height actions
        WSCEst_Rj       = initWSC*ones(1,nRj);          %   Action value estimation vector for orbit radius actions
        WSCN_Angle      = zeros(1,nAng);                %   Vector to store angle action ocurrences
        WSCN_H          = zeros(1,nH);                  %   Vector to store height action ocurrences
        WSCN_Rj         = zeros(1,nRj);                 %   Vector to store orbit radius action ocurrences

        WSC_UAV = 0;                                    %   WSC initial value
        WSC_Ang = angleUAV(end);               %   Angle random initialization
        WSC_H   = hUAV(end);                   %   Height random initialization
        WSC_Rj  = rjUAV(end);                  %   Orbit Radius random initialization

        tic
        for i=1:nLoops
            if mod(i,nRLB+1)==0
                % This is a positioning block
                % No RL is performed, and Performance variables are taken at
                % this block only
                t1 = toc;                                   % Time from previous SIB end to this Positioning block
                iSIB = i/(nRLB+1);                          % Number of SIB

                % Obtain new position indeces
                [indAng, indH, indRj] = getNewPos(WSCEst_Angle, WSCEst_H, WSCEst_Rj);   

                % UAV Movement
                WSC_Rj_Old = WSC_Rj;                                                    
                dH      = abs( WSC_H - hUAV(indH) );                                    % Height movement distance
                dRj     = abs( WSC_Rj - rjUAV(indRj) );                                 % Orbit radius movement distance
                dAng    = min(WSC_Rj,WSC_Rj_Old)*abs( WSC_Ang - angleUAV(indAng) )/2;	% Angle movement distance (through lesser orbit)
                dtJ     = (dAng + dH + dRj)/v_J;                                        % Flight time

                % Update UAV positioning values
                WSC_Ang = angleUAV(indAng);        
                WSC_H   = hUAV(indH);
                WSC_Rj  = rjUAV(indRj);

                % Update UAV position
                UAVs = setNewPos_N(nUAV, WSC_Ang, WSC_H, WSC_Rj);

                % Compute secrecy performance metric
                S_SIB(iRLB,iSIB) = S_SIB(iRLB,iSIB) +  dtRL*nRLB*3*WSC_UAV/nMC;

                % Compute energy metric
                E_SIB(iRLB,iSIB) = E_SIB(iRLB,iSIB) +  (E_RX + E_ACK + dtJ*P_Mov)/nMC;

                % Update UAV WSC value
                [WSC_UAV, DeltaComp] = computeWSC_NUAV(A, E, UAVs, dAB_R, gammaA, gammaJ, channelParam, boolFade );
                WSC_UAV = WSC_UAV/(nTheta*nR);

                % Store control values per MC loop
                MC_WSC(iRLB,iMC,iSIB)    = WSC_UAV;
                MC_Ang(iRLB,iMC,iSIB)    = WSC_Ang;               %   Stored Angles per SIB for all MC loops
                MC_H(iRLB,iMC,iSIB)      = WSC_H;                 %   Stored Heights per SIB for all MC loops
                MC_Rj(iRLB,iMC,iSIB)     = WSC_Rj;                %   Stored Orbit Radius per SIB for all MC loops

                fprintf('K: %.1f \t\t MCLoop: %i of %i \t\t Positioning Block: %i of %i\t\tTime: %.3f\t\t WSC: %.3f\n',KbVal,iMC,nMC,iSIB,nSIB, t1,WSC_UAV);
                tic
                continue;
            end
            % If it is not a Positioning Block, it's an RL block
            % A single RL iteration is performed by each RL variable in
            % succession (angle, height, orbit radius) for a single dAB
            % estimate

            % dAB estimation and parameter computation
            dAB = normrnd(dAB_R,sigAB);                             %   CSI estimate

            %   RL iteration for angle (greedy on h and Rj)
            [WSCEst_Angle, WSCN_Angle]  = computeRLBlock_UCB_N(WSCEst_Angle, WSCEst_H, WSCEst_Rj, WSCN_Angle, angleUAV, hUAV, rjUAV, ...
                                                A, E, dAB, gammaA, gammaJ, c, channelParam, i, 1, alpha, nUAV, boolFade);

            %   RL iteration for H (greedy on angle and Rj)
            [WSCEst_H, WSCN_H]          = computeRLBlock_UCB_N(WSCEst_Angle, WSCEst_H, WSCEst_Rj, WSCN_H, angleUAV, hUAV, rjUAV, ...
                                                A, E, dAB, gammaA, gammaJ, c, channelParam, i, 2, alpha, nUAV, boolFade);

            %   RL iteration for Rj (greedy on angle and H)
            [WSCEst_Rj, WSCN_Rj]        = computeRLBlock_UCB_N(WSCEst_Angle, WSCEst_H, WSCEst_Rj, WSCN_Rj, angleUAV, hUAV, rjUAV, ...
                                                A, E, dAB, gammaA, gammaJ, c, channelParam, i, 3, alpha, nUAV, boolFade);
        end
    end
    save(['SIB_EX_nK_GenConst-v0F_', num2str(iRLB), '.mat'])
end

save('SIB_EX_nK_nUAV2-v10-fade.mat')
%% Multi Varying nUAV
clc, clear, close all

% Environment parameters (Urban)

phi = 9.61;         % Environmental constant
omega = 0.16 ;      % Environmental constant
alpha = 0.3;        % Ground Path Loss exponent
alpha_AG = 0.3;     % Air-to-Ground Path Loss Exponent
ne_LOS = 1.0;       % Air-to-Ground LOS attenuation
ne_NLOS = 20;       % Air-to-Ground NLOS attenuation

KbVal = 5;
KeVal = 5;

% Channel parameters

m = 3;              % Number of parallel channels
sigma = 1/sqrt(2);  % Noise std dev of each component (real and imag) of every parallel channel

Rs = 1;

channelParam =  [   phi,...
                    omega,...
                    alpha,...
                    alpha_AG,...
                    ne_LOS,...
                    ne_NLOS,...
                    Rs,...
                    KbVal,...
                    KeVal];


% Positions ***************************************************************

    % Alice
nA = 1;             %   Number of Alices
A = [0,0,0];        %   Position of Alice (zero point)
gammaA = 200;        %   Alice Tx SNR

% -------------------------------------------------------------------------

    % Bob
nB = 1;             %   Number of Bobs
dAB_R = 20;        	%   Real distance between A and B
sigAB = dAB_R/10;
B = [dAB_R,0,0];      %   Position of Bob

% -------------------------------------------------------------------------

    % Eve
nR = 100;           %   Number of radial points
nTheta = 360;       %   Number of angular points
nE = nR*nTheta;     %   Number of Eves

rLow = 0.1;         %   Lowest radius of Eve
rHigh = 50;         %   Highest radius of Eve
thetaLow = 0;       %   Lowest angle of Eve
thetaHigh = 2*pi;   %   Highest radius of Eve

rangeR = linspace(rLow,rHigh,nR);                       % Points in Radial dimension
rangeTheta = linspace(thetaLow,thetaHigh,nTheta);       % Points in Angular dimension

thetat  =  repmat(rangeTheta,1,nR);
rt  = (repmat(rangeR',1,nTheta)).';

E = [rt(:).*cos(thetat(:)), rt(:).*sin(thetat(:)) , zeros(nR*nTheta,1)];             % Eves' position (rectangle coordinates)

% -------------------------------------------------------------------------

    % UAV
Nu = 2;                             %   Number of simultaneous UAVs
nUAV_v = [2 4 8];
NnUAV = length(nUAV_v);

nAng = 90;                          %   Angle discretization level (opening angle)  -> Number of Angle Actions
nH = 30;                            %   Height discretization level for both UAVs   -> Number of Height Actions
nRj = 30;                           %   Orbit radius discretization level           -> Number of radius Actions

hUAV        = linspace(100,200,nH);                    %   Possible height actions
rjUAV       = linspace(1,rHigh,nRj);                %   Possible orbit radius actions

% -------------------------------------------------------------------------

% *************************************************************************


%   UAV technical parameters [dummy values, research due]

E_RX = 0.1;                 % Energy it takes for UAV to receive positional information [dummy value].
E_ACK = 0.1;                % Energy it takes for UAV to send ACK to Alice [dummy value].
P_Mov = 80;                 % Power it takes UAV to maneuver [dummy value].
v_J = 2;                    % Speed at which UAV moves [dummy value].

%   SIB k-Armed Bandits
nSIB = 20;                  % Number of SIBs to be considered
nRLB = 10; 
nLoops = nSIB*(nRLB+1);     % Number of loops for action choosing (positioning every nRLB+1 loops)

nRL = 1;                    % Number of RL steps per estimate
epsilon = 0;                % Epsilon for epsilon-greedy choice
initWSC = 0;                % Optimistic initial action values
c = 0.3;                    % Exploration parameter for UCB
alpha = 0.5;                  % Step size (0: uniform average)

%   Performance Variables

nMC = 1e2;                              %   Number of Monte Carlo loops for performance evaluation
E_SIB = zeros(NnUAV,nSIB);                  %   Energy used in a given SIB
S_SIB = zeros(NnUAV,nSIB);                  %   Cummulative secrecy over a given SIB

MC_WSC = zeros(NnUAV,nMC,nSIB);               %   Stored WSC per SIB for all MC loops
MC_Ang = zeros(NnUAV,nMC,nSIB);               %   Stored Angles per SIB for all MC loops
MC_H = zeros(NnUAV,nMC,nSIB);                 %   Stored Heights per SIB for all MC loops
MC_Rj = zeros(NnUAV,nMC,nSIB);                %   Stored Orbit Radius per SIB for all MC loops

dtRL = 0.1*3;                             % DUMMY TIME, to get consistent results, not depending on the processor or machine
                                        % Time per RL iteration (the smallest RL iteration)

boolFade = 1;               % Use fading


for iRLB = 1:NnUAV
    nUAV = NnUAV(iRLB);
    gammaJ = 0.1*gammaA/nUAV;                 %   UAVs Jamming SNR
    angleUAV    = linspace(0,2*pi/(nUAV-1),nAng);       %   Possible angle actions (opening angles)

    % Exhaustive Search
    tic
    fprintf('Exhaustive search!\n')
    WSC_ES(iRLB,:,:,:) = exhaustiveSearchWSC(A, E, nUAV, angleUAV, hUAV, rjUAV, dAB_R, gammaA, gammaJ, channelParam, boolFade );
    fprintf('\n\n Exhaustive search ended!')
    toc
    for iMC = 1:nMC
        % Initialization
        WSCEst_Angle    = initWSC*ones(1,nAng);         %   Action value estimation vector for angle actions
        WSCEst_H        = initWSC*ones(1,nH);           %   Action value estimation vector for height actions
        WSCEst_Rj       = initWSC*ones(1,nRj);          %   Action value estimation vector for orbit radius actions
        WSCN_Angle      = zeros(1,nAng);                %   Vector to store angle action ocurrences
        WSCN_H          = zeros(1,nH);                  %   Vector to store height action ocurrences
        WSCN_Rj         = zeros(1,nRj);                 %   Vector to store orbit radius action ocurrences

        WSC_UAV = 0;                                    %   WSC initial value
        
        WSC_Ang = angleUAV(end);               %   Angle random initialization
        WSC_H   = hUAV(end);                   %   Height random initialization
        WSC_Rj  = rjUAV(end);                  %   Orbit Radius random initialization

        tic
        for i=1:nLoops
            if mod(i,nRLB+1)==0
                % This is a positioning block
                % No RL is performed, and Performance variables are taken at
                % this block only
                t1 = toc;                                   % Time from previous SIB end to this Positioning block
                iSIB = i/(nRLB+1);                          % Number of SIB

                % Obtain new position indeces
                [indAng, indH, indRj] = getNewPos(WSCEst_Angle, WSCEst_H, WSCEst_Rj);   

                % UAV Movement
                WSC_Rj_Old = WSC_Rj;                                                    
                dH      = abs( WSC_H - hUAV(indH) );                                    % Height movement distance
                dRj     = abs( WSC_Rj - rjUAV(indRj) );                                 % Orbit radius movement distance
                dAng    = min(WSC_Rj,WSC_Rj_Old)*abs( WSC_Ang - angleUAV(indAng) )/2;	% Angle movement distance (through lesser orbit)
                dtJ     = (dAng + dH + dRj)/v_J;                                        % Flight time

                % Update UAV positioning values
                WSC_Ang = angleUAV(indAng);        
                WSC_H   = hUAV(indH);
                WSC_Rj  = rjUAV(indRj);

                % Update UAV position
                UAVs = setNewPos_N(nUAV, WSC_Ang, WSC_H, WSC_Rj);

                % Compute secrecy performance metric
                S_SIB(iRLB,iSIB) = S_SIB(iRLB,iSIB) +  dtRL*nRLB*3*WSC_UAV/nMC;

                % Compute energy metric
                E_SIB(iRLB,iSIB) = E_SIB(iRLB,iSIB) +  (E_RX + E_ACK + dtJ*P_Mov)/nMC;

                % Update UAV WSC value
                [WSC_UAV, DeltaComp] = computeWSC_NUAV(A, E, UAVs, dAB_R, gammaA, gammaJ, channelParam, boolFade );
                WSC_UAV = WSC_UAV/(nTheta*nR);

                % Store control values per MC loop
                MC_WSC(iRLB,iMC,iSIB)    = WSC_UAV;
                MC_Ang(iRLB,iMC,iSIB)    = WSC_Ang;               %   Stored Angles per SIB for all MC loops
                MC_H(iRLB,iMC,iSIB)      = WSC_H;                 %   Stored Heights per SIB for all MC loops
                MC_Rj(iRLB,iMC,iSIB)     = WSC_Rj;                %   Stored Orbit Radius per SIB for all MC loops

                fprintf('nUAV: %.1f \t\t MCLoop: %i of %i \t\t Positioning Block: %i of %i\t\tTime: %.3f\t\t WSC: %.3f\n',nUAV,iMC,nMC,iSIB,nSIB, t1,WSC_UAV);
                tic
                continue;
            end
            % If it is not a Positioning Block, it's an RL block
            % A single RL iteration is performed by each RL variable in
            % succession (angle, height, orbit radius) for a single dAB
            % estimate

            % dAB estimation and parameter computation
            dAB = normrnd(dAB_R,sigAB);                             %   CSI estimate

            %   RL iteration for angle (greedy on h and Rj)
            [WSCEst_Angle, WSCN_Angle]  = computeRLBlock_UCB_N(WSCEst_Angle, WSCEst_H, WSCEst_Rj, WSCN_Angle, angleUAV, hUAV, rjUAV, ...
                                                A, E, dAB, gammaA, gammaJ, c, channelParam, i, 1, alpha, nUAV, boolFade);

            %   RL iteration for H (greedy on angle and Rj)
            [WSCEst_H, WSCN_H]          = computeRLBlock_UCB_N(WSCEst_Angle, WSCEst_H, WSCEst_Rj, WSCN_H, angleUAV, hUAV, rjUAV, ...
                                                A, E, dAB, gammaA, gammaJ, c, channelParam, i, 2, alpha, nUAV, boolFade);

            %   RL iteration for Rj (greedy on angle and H)
            [WSCEst_Rj, WSCN_Rj]        = computeRLBlock_UCB_N(WSCEst_Angle, WSCEst_H, WSCEst_Rj, WSCN_Rj, angleUAV, hUAV, rjUAV, ...
                                                A, E, dAB, gammaA, gammaJ, c, channelParam, i, 3, alpha, nUAV, boolFade);
        end
    end
    save(['SIB_EX_nUAVN_GenConst-v0F_', num2str(iRLB), '.mat'])
end

save('SIB_EX_nUAVN-v10-fade.mat')

%% Multi Varying c
clc, clear, close all

% Environment parameters (Urban)

phi = 9.61;         % Environmental constant
omega = 0.16 ;      % Environmental constant
alpha = 0.3;        % Ground Path Loss exponent
alpha_AG = 0.3;     % Air-to-Ground Path Loss Exponent
ne_LOS = 1.0;       % Air-to-Ground LOS attenuation
ne_NLOS = 20;       % Air-to-Ground NLOS attenuation
KbVal = 5;          % K parameter for Bob's channels
KeVal = 5;          % K parameter for Eve's channels


% Channel parameters

m = 3;              % Number of parallel channels
sigma = 1/sqrt(2);  % Noise std dev of each component (real and imag) of every parallel channel

Rs = 1;


channelParam =  [   phi,...
                    omega,...
                    alpha,...
                    alpha_AG,...
                    ne_LOS,...
                    ne_NLOS,...
                    Rs,...
                    KbVal,...
                    KeVal];

% Positions ***************************************************************

    % Alice
nA = 1;             %   Number of Alices
A = [0,0,0];        %   Position of Alice (zero point)
gammaA = 200;        %   Alice Tx SNR

% -------------------------------------------------------------------------

    % Bob
nB = 1;             %   Number of Bobs
dAB_R = 20;        	%   Real distance between A and B
sigAB = dAB_R/10;          %   Unreliability of B's position
B = [dAB_R,0,0];      %   Position of Bob

% -------------------------------------------------------------------------

    % Eve
nR = 100;           %   Number of radial points
nTheta = 360;       %   Number of angular points
nE = nR*nTheta;     %   Number of Eves

rLow = 0.1;         %   Lowest radius of Eve
rHigh = 50;         %   Highest radius of Eve
thetaLow = 0;       %   Lowest angle of Eve
thetaHigh = 2*pi;   %   Highest radius of Eve

rangeR = linspace(rLow,rHigh,nR);                       % Points in Radial dimension
rangeTheta = linspace(thetaLow,thetaHigh,nTheta);       % Points in Angular dimension

thetat  =  repmat(rangeTheta,1,nR);
rt  = (repmat(rangeR',1,nTheta)).';

E = [rt(:).*cos(thetat(:)), rt(:).*sin(thetat(:)) , zeros(nR*nTheta,1)];             % Eves' position (rectangle coordinates)

% -------------------------------------------------------------------------

    % UAV
Nu = 2;                             %   Number of simultaneous UAVs
nUAV = 2;

nAng = 90;                          %   Angle discretization level (opening angle)  -> Number of Angle Actions
nH = 30;                            %   Height discretization level for both UAVs   -> Number of Height Actions
nRj = 30;                           %   Orbit radius discretization level           -> Number of radius Actions

gammaJ = 0.1*gammaA/nUAV;                 %   UAVs Jamming SNR

angleUAV    = linspace(0,2*pi/(nUAV-1),nAng);       %   Possible angle actions (opening angles)
hUAV        = linspace(100,200,nH);                    %   Possible height actions
rjUAV       = linspace(1,rHigh,nRj);                %   Possible orbit radius actions

% -------------------------------------------------------------------------

% *************************************************************************


%   UAV technical parameters [dummy values, research due]

E_RX = 0.1;                 % Energy it takes for UAV to receive positional information [dummy value].
E_ACK = 0.1;                % Energy it takes for UAV to send ACK to Alice [dummy value].
P_Mov = 80;                 % Power it takes UAV to maneuver [dummy value].
v_J = 2;                    % Speed at which UAV moves [dummy value].

%   SIB k-Armed Bandits
nSIB = 20;                  % Number of SIBs to be considered
nRLB = 10; 
nLoops = nSIB*(nRLB+1);     % Number of loops for action choosing (positioning every nRLB+1 loops)

nRL = 1;                    % Number of RL steps per estimate
epsilon = 0;                % Epsilon for epsilon-greedy choice
initWSC = 0;                % Optimistic initial action values

cV = [0.01 0.1 0.3];                  % Exploration parameter for UCB
nC = length(cV);

%c = 0.3;                    % Exploration parameter for UCB
alpha = 0.5;                  % Step size (0: uniform average)

%   Performance Variables

nMC = 1e2;                              %   Number of Monte Carlo loops for performance evaluation
E_SIB = zeros(nC,nSIB);                  %   Energy used in a given SIB
S_SIB = zeros(nC,nSIB);                  %   Cummulative secrecy over a given SIB

MC_WSC = zeros(nC,nMC,nSIB);               %   Stored WSC per SIB for all MC loops
MC_Ang = zeros(nC,nMC,nSIB);               %   Stored Angles per SIB for all MC loops
MC_H = zeros(nC,nMC,nSIB);                 %   Stored Heights per SIB for all MC loops
MC_Rj = zeros(nC,nMC,nSIB);                %   Stored Orbit Radius per SIB for all MC loops

dtRL = 0.1*3;                             % DUMMY TIME, to get consistent results, not depending on the processor or machine
                                        % Time per RL iteration (the smallest RL iteration)

boolFade = 1;               % Use fading

% Exhaustive Search
tic
fprintf('Exhaustive search!\n')
WSC_ES = exhaustiveSearchWSC(A, E, nUAV, angleUAV, hUAV, rjUAV, dAB_R, gammaA, gammaJ, channelParam, boolFade );
fprintf('\n\n Exhaustive search ended!')
toc



for iRLB = 1:nC
    c = cV(iRLB);
    for iMC = 1:nMC
        % Initialization
        WSCEst_Angle    = initWSC*ones(1,nAng);         %   Action value estimation vector for angle actions
        WSCEst_H        = initWSC*ones(1,nH);           %   Action value estimation vector for height actions
        WSCEst_Rj       = initWSC*ones(1,nRj);          %   Action value estimation vector for orbit radius actions
        WSCN_Angle      = zeros(1,nAng);                %   Vector to store angle action ocurrences
        WSCN_H          = zeros(1,nH);                  %   Vector to store height action ocurrences
        WSCN_Rj         = zeros(1,nRj);                 %   Vector to store orbit radius action ocurrences

        WSC_UAV = 0;                                    %   WSC initial value

        WSC_Ang = angleUAV(end);               %   Angle random initialization
        WSC_H   = hUAV(end);                   %   Height random initialization
        WSC_Rj  = rjUAV(end);                  %   Orbit Radius random initialization

        tic
        for i=1:nLoops
            if mod(i,nRLB+1)==0
                % This is a positioning block
                % No RL is performed, and Performance variables are taken at
                % this block only
                t1 = toc;                                   % Time from previous SIB end to this Positioning block
                iSIB = i/(nRLB+1);                          % Number of SIB

                % Obtain new position indeces
                [indAng, indH, indRj] = getNewPos(WSCEst_Angle, WSCEst_H, WSCEst_Rj);   

                % UAV Movement
                WSC_Rj_Old = WSC_Rj;                                                    
                dH      = abs( WSC_H - hUAV(indH) );                                    % Height movement distance
                dRj     = abs( WSC_Rj - rjUAV(indRj) );                                 % Orbit radius movement distance
                dAng    = min(WSC_Rj,WSC_Rj_Old)*abs( WSC_Ang - angleUAV(indAng) )/2;	% Angle movement distance (through lesser orbit)
                dtJ     = (dAng + dH + dRj)/v_J;                                        % Flight time

                % Update UAV positioning values
                WSC_Ang = angleUAV(indAng);        
                WSC_H   = hUAV(indH);
                WSC_Rj  = rjUAV(indRj);

                % Update UAV position
                UAVs = setNewPos_N(nUAV, WSC_Ang, WSC_H, WSC_Rj);

                % Compute secrecy performance metric
                S_SIB(iRLB,iSIB) = S_SIB(iRLB,iSIB) +  dtRL*nRLB*3*WSC_UAV/nMC;

                % Compute energy metric
                E_SIB(iRLB,iSIB) = E_SIB(iRLB,iSIB) +  (E_RX + E_ACK + dtJ*P_Mov)/nMC;

                % Update UAV WSC value
                [WSC_UAV, DeltaComp] = computeWSC_NUAV(A, E, UAVs, dAB_R, gammaA, gammaJ, channelParam, boolFade );
                WSC_UAV = WSC_UAV/(nTheta*nR);

                % Store control values per MC loop
                MC_WSC(iRLB,iMC,iSIB)    = WSC_UAV;
                MC_Ang(iRLB,iMC,iSIB)    = WSC_Ang;               %   Stored Angles per SIB for all MC loops
                MC_H(iRLB,iMC,iSIB)      = WSC_H;                 %   Stored Heights per SIB for all MC loops
                MC_Rj(iRLB,iMC,iSIB)     = WSC_Rj;                %   Stored Orbit Radius per SIB for all MC loops

                fprintf('c: %.3f \t\t MCLoop: %i of %i \t\t Positioning Block: %i of %i\t\tTime: %.3f\t\t WSC: %.3f\n',c,iMC,nMC,iSIB,nSIB, t1,WSC_UAV);
                tic
                continue;
            end
            % If it is not a Positioning Block, it's an RL block
            % A single RL iteration is performed by each RL variable in
            % succession (angle, height, orbit radius) for a single dAB
            % estimate

            % dAB estimation and parameter computation
            dAB = normrnd(dAB_R,sigAB);                             %   CSI estimate

            %   RL iteration for angle (greedy on h and Rj)
            [WSCEst_Angle, WSCN_Angle]  = computeRLBlock_UCB_N(WSCEst_Angle, WSCEst_H, WSCEst_Rj, WSCN_Angle, angleUAV, hUAV, rjUAV, ...
                                                A, E, dAB, gammaA, gammaJ, c, channelParam, i, 1, alpha, nUAV, boolFade);

            %   RL iteration for H (greedy on angle and Rj)
            [WSCEst_H, WSCN_H]          = computeRLBlock_UCB_N(WSCEst_Angle, WSCEst_H, WSCEst_Rj, WSCN_H, angleUAV, hUAV, rjUAV, ...
                                                A, E, dAB, gammaA, gammaJ, c, channelParam, i, 2, alpha, nUAV, boolFade);

            %   RL iteration for Rj (greedy on angle and H)
            [WSCEst_Rj, WSCN_Rj]        = computeRLBlock_UCB_N(WSCEst_Angle, WSCEst_H, WSCEst_Rj, WSCN_Rj, angleUAV, hUAV, rjUAV, ...
                                                A, E, dAB, gammaA, gammaJ, c, channelParam, i, 3, alpha, nUAV, boolFade);
        end
    end
    save(['SIB_EX_nC_GenConst-v0F_', num2str(iRLB), '.mat'])
end

save('SIB_EX_nC_nUAV2-v10-fade.mat')

🎉3 参考文献 

文章中一些内容引自网络,会注明出处或引用为参考文献,难免有未尽之处,如有不妥,请随时联系删除。(文章内容仅供参考,具体效果以运行结果为准)

🌈Matlab代码实现

资料获取,更多粉丝福利,MATLAB|Simulink|Python资源获取

                                                           在这里插入图片描述

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值