基于连续相位负载调制的单输入宽带混合Doherty功率放大器设计-从理论到ADS版图
最近开始搞大论文了,Doherty方面对最新的一些的技术看的比较少,找几个牛逼的学习一下下,虽然最后可能也用不上。已经完成了理论的推导和ADS版图仿真,在全部介绍完成后会一起分享出来。
原文: Single-Input Broadband Hybrid Doherty Power Amplifiers Design Relying on a Phase Sliding-Mode of the Load Modulation Scheme
发表于APRIL 2023,在微波顶刊IEEE T MTT上面,使用的GAN CGH40010F
全部工程下载:https://download.csdn.net/download/weixin_44584198/88795215
下载前阅读文章4、文件解释(作者我复现了一周,而且是MTT的最新理论,当个硕士毕业论文的一大章戳戳有余嘞;所以跟感兴趣的金主大大们要点饭钱,穷学生私信我,可以打个折)
当然也可以自己复现,毕竟论文就在那边嘞
目录
- 基于连续相位负载调制的单输入宽带混合Doherty功率放大器设计-从理论到ADS版图
- 0、实现效果展示
- 0.1、版图展示
- 0.2、单音信号版图仿真结果
- 0.3、调制信号与DPD版图仿真结果
- 1、文章的核心思路
- 1.1、提要
- 1.2、核心思路
- 2、推导与图表复现
- 2.1 Design of the Phase Sliding-Mode Output Combiner
- 2.2 输入移相器设计
- 2.3 理想的性能表现
- 3、ADS工程构建与测试
- 3.1、功分器
- 3.2、输入匹配设计
- 3.3、输出合路器设计
- 3.4、后匹配电路设计
- 3.5、联合仿真
- 4、文件解释(下载资源前必看)
- 4.1 环境
- 4.2 文件解释
- 5、Matlab代码
0、实现效果展示
原来的作者在论文中贴出了微带电路图和原理图,但是原理图和微带电路图实际对不上,而且基于作者给出的原理图进行仿真结果对不上,因此我使用作者提出的理论自己设计了一份,电路结构有所不同但是效果基本类似的。下面所提及的电路、效果,如非特别提及,都是我使用作者理论自己设计并在ADS中版图仿真的结果。
0.1、版图展示
整体的版图如下(偏置边上的大块接地铜皮是给焊接电容用的,使用的板材和作者文章中的一致):
0.2、单音信号版图仿真结果
基于上面的仿真版图,进行谐波调谐的仿真(仿真使用的是CRCW0603电阻和村田的GRM18电容在ADS中的模型),运行HB1ToneGComp1swp_Doherty原理图:
可以看到,在2.4-3.2GHz的范围内,饱和功率在43-44dBm之间;饱和效率在63.3%-68%之间;如果认为饱和功率为44dBm(此处近似一下,不然解释有点烦),回退6dB的效率为52.5%-56.7%,回退8dB的效率为45.4%-49%。
在2.2-3.4GHz的范围内,回退6dB的效率在45%以上,回退8dB的效率在41%以上,饱和效率在61%以上(图片中就不展示了,不是设计目标)。
0.3、调制信号与DPD版图仿真结果
运行DPD_Test或者直接打开dds图,ADS可能无法直接运行,详见番外10:使用ADS对射频功率放大器进行非线性测试2(使用带宽20MHz的64QAM信号进行ACLR、EVM、CCDF测试)。此处直接放出在2.4GHz的使用64QAM调制的20MHz带宽的LTE信号的结果:
由图可见,经过DPD后的ACLR小于-45dBc,符合相关通信要求,DPD对链路的改善为17dBc。使用调制信号测试的平均漏极效率为52.5%以上,平均PAE为49.5%左右。其余详细也见上图。
1、文章的核心思路
1.1、提要
理想架构的Doherty功率放大器理论与仿真中介绍了传统DPA的回退思路,理想架构的非对称高回退Doherty功率放大器理论与仿真中也介绍了高回退的非对称的DPA的设计理论,但是还未对其宽带性能进行讨论。
先看看传统DPA的限制,DPA进行负载调制的基础就是四分之一波长线的阻抗变换作用,但众所周知,四分之一波长线带宽很窄,只在单一频点上效果较好(这一点可以使用理想架构的Doherty功率放大器理论与仿真中的模板进行仿真,我为了精简就不展开了):
为了解决这个问题,作者提出了基于连续相位负载调制的单输入宽带混合Doherty功率放大器,其思路来源于Doherty–Chireix(异相)连续模式,那篇文章我没看,好像也是Trans上的,但是我暂时没发现这个理论和Outphasing的类似之处。
1.2、核心思路
作者提出了一种新的结构,其和上面传统的结构的区别就是峰值(辅助)功放的输出也接了一个微带线。
更为具体的是,对于传统DPA分析,四分之一波长线的相位延迟是固定的,对于所设计的频点是90°,对于其他频率点也会有固定成规律的延迟,因此实际上 θ \theta θ是关于频率的变量。
对于非中心频率,四分之一波长线所提供的 θ \theta θ不是90°,无法满足负载调制关系(下面微带线是2.8GHz的四分之一波长线的延迟):
核心的问题就是四分之一波长线所提供的 θ \theta θ是频率的变量,如果能够找到一个结构能够在宽频带内实现同样的功能,那么宽带DPA的设计就不是难题了,当然找不到。
本文的贡献就是,提出了一种新的DPA设计方法,需要将峰值功放、载波功放后的微带结构设计成和频率相关的特定相位延迟,这样就能在宽带内实现DPA结构。在这种理论中,负载RL的值不是固定的,也是和频率变化的函数。
以文章中给出的结果为例,观察绿色方框,分别给出理论的微带延迟和频率的关系,其中 θ a \theta_a θa是辅助功放(Auxiliary)的延迟表达, θ m \theta_m θm是主功放(Main)的微带结构应该具备的延迟;负载所需阻抗大概集中在10欧姆附近,这样随频率变化的负载阻抗设计并不难,因为我们有后匹配网络,设计成啥样都行:
此外, θ \theta θ是主、辅功放的相位延迟,会在输入栅极使用微带线进行补偿。主功放的微带线特性阻抗Zm被设置在24.4欧姆(Ropt最佳B类阻抗,对于CGH40010F这个管子选择24、30、36我都见过),而辅助功放接的微带线特性阻抗的实部Rap是一个随频率变化的变量。
2、推导与图表复现
这部分仅仅对想要深入了解论文的人有用。基于文章中的公式推导,我编写了Matlab代码并对理论进行了仿真ADS绘图等等。
2.1 Design of the Phase Sliding-Mode Output Combiner
作者首先使用了定义:辅助 PA 的峰值与回退漏极电压之比(参考论文Novel outphasing power amplifiers designed with an analytic generalized Doherty–Chireix continuum theory、Wideband two-way hybrid Doherty outphasing power amplifier):
K v a = ∣ V a p ∣ / ∣ V a b ∣ K_{va}=|V_{ap}|/|V_{ab}| Kva=∣Vap∣/∣Vab∣
值得强调一点,辅助功放在回退点不工作,其实质是有电压没电流的,不清楚的可以看看DPA理想仿真里面的介绍(理想架构的Doherty功率放大器理论与仿真),这边贴图(左边电流图,右边电压图):
由此对于经典的DPA来说, K v a = ∣ V a p ∣ / ∣ V a b ∣ = 2 K_{va}=|V_{ap}|/|V_{ab}|=2 Kva=∣Vap∣/∣Vab∣=2。因为回退时峰值功放不工作,没有电流,因此辅助 PA 的峰值与回退漏极电流之比是无穷:
K i a = ∣ I a p ∣ / ∣ I a b ∣ = ∞ K_{ia}=|I_{ap}|/|I_{ab}|=\infty Kia=∣Iap∣/∣Iab∣=∞
结合一些其他论文的式子,可以得出功率分配比的计算公式(功率分配比在DPA的概念参考:理想架构的非对称高回退Doherty功率放大器理论与仿真):
n ( ω ) = K v a ( ω ) O B O − 1 O B O + K v a ( ω ) n(\omega)=K_{va}(\omega)\frac{\mathrm{OBO}-1}{\mathrm{OBO}+K_{va}(\omega)} n(ω)=Kva(ω)OBO+Kva(ω)OBO−1
值得注意的是OBO代表的回退比,对于传统的DPA是6dB,但是计算的时候需要把它单位转换为正常的不用对数,如:
O B O = 1 0 ( O B O d B / 10 ) OBO=10^{(OBO_{dB}/10)} OBO=10(OBOdB/10)
对于传统的DPA是6dB,那么OBO就是10^0.6=4。那么传统DPA的功率分配比是:
n ( ω ) = K v a ( ω ) O B O − 1 O B O + K v a ( ω ) = 2 ∗ ( 4 − 1 ) / ( 4 + 2 ) = 1 n(\omega)=K_{va}(\omega)\frac{\mathrm{OBO}-1}{\mathrm{OBO}+K_{va}(\omega)}=2*(4-1)/(4+2)=1 n(ω)=Kva(ω)OBO+Kva(ω)OBO−1=2∗(4−1)/(4+2)=1
就是1:1等分的DPA架构啦。
其余的就是给出公式和解方程,在自己Matlab中进行并画出和图3一样的效果,最终的结果如下,代码在文章后面放出来(结果和作者原版有细小的差别,猜测是作者给出的Kva数据是近似过的):
这个图实际上就是作者计算出来的各个微带结构要具备多少的相位。从图中也可以看到,一共有四种选择,作者接下来就是花篇幅去解释如何选择合理的,对于θm和θa,应选择两个正号(θm1和θa1)或两个负号(θm4和θa4)来减少两个TL(θm和θa)的物理长度。
2.2 输入移相器设计
上面还有两个重要的角度参数 θ p \theta_{p} θp和 θ b \theta_{b} θb,其中 θ b \theta_{b} θb是回退时输入所需要具备的异相角, θ p \theta_{p} θp是峰值时输入所需要具备的异相角。
但是,对于双输入的Outphasing架构,异相角与功率相关并不难设计。然而,为了实现单输入 PA,如果要使用无源输入移相器电路来实现,则异相角必定与功率无关。对此,作者使用了近似把这两个角度设置为相等。这样的话一个有8个数值可以选择,也就是 θ p 1 − θ p 4 \theta_{p1}-\theta_{p4} θp1−θp4, θ b 1 − θ b 4 \theta_{b1}-\theta_{b4} θb1−θb4。
对于这8个,也是要选择小的节省微带线长度,也就是θb1、θb3、θp2、θp4这四个里面选(参考上面的图)。但是在2.1中,我们选择同号的来减少θm和θa的物理长度,就是使用的1、4,因此我们异相角的选择只有θb1和θp4了。
作者接下来的推导我没有看懂,从图中看θb1和θp4的数值是相同的,那么选择任何一个应该都没有区别啊,看懂的小伙伴可以评论区告诉我。我直接展示一下作者的结果了(我自己代码运行得到的两个都是紫色那个线条):
作者:可以看到,选择θp4可以获得更加稳定的饱和输出功率,效果更好;作者也进行了ADS理论仿真,结果和上图基本一致:
2.3 理想的性能表现
作者提出了在不同频率点的多种角度、阻抗组合,可以在宽带内实现DPA架构,下面使用理想的ADS电流源仿真观察其理想性能,参数表如下图:
为了简介,我们只对Kva=1.75的情况在ADS中进行仿真,其余的情况都是类似的,构建ADS电路图(图中器件参数都是我代码计算出来的,和作者结果相差不大,这个理想仿真工程文件包含在最上方链接):
理想的仿真结果如下所示,由此可见作者提出的理论是可行的,能够实现8dB范围的回退。
但是由于OBO固定,每个Kva数值下的功率分配比不是固定的(上图大概是1:1.15的分配比的样子),分配比计算公式为: n ( ω ) = K v a ( ω ) O B O − 1 O B O + K v a ( ω ) n(\omega)=K_{va}(\omega)\frac{\mathrm{OBO}-1}{\mathrm{OBO}+K_{va}(\omega)} n(ω)=Kva(ω)OBO+Kva(ω)OBO−1
作者Fig6原图如下,展示了各个Kva值的回退性能,可以从右图看到似乎Kva越大饱和输出功率越高,这是因为Kva越大分配越不平衡,在载波功放饱和时峰值功放过饱和了,因此输出功率大。可以看到Kva=1.5分配比约为1:1,输出功率为2;Kva=2.5分配比约为1:1.5,总输出功率为2.5,相差-10*log(2/2.5)=1dB,和图中一致:
3、ADS工程构建与测试
3.1、功分器
作者图里面的功分器没有具体结构,需要自己设计,使用ADS自带的微带线功分器设计工具,设计结果如下(特意往高频率设计了点,差别不大):
3.2、输入匹配设计
原来作者提供的输入匹配的原理图和其实物图对不上,只能自己设计了,设计结构如下(原始数据来自源牵引,源牵引教程参考ADS功放设计之负载牵引与源牵引):
3.3、输出合路器设计
依据作者提出的理论,输出合路器的结构需要满足一定的要求,也就是峰值功放和载波功放输出微带结构的延迟需要和理论推导一致,作者给出其理论和实际电路图实现的对比:
使用作者给出的电路原理图,结果是基本一致的(但是原理图和版图不一致),首先看看主功放后面电路实现的结果 :
此处我表述错误,是我没有考虑到封装参数,那就是作者给出的原理图有误-=-
3.4、后匹配电路设计
对于不同频段,要求的负载RL不一致,但是总体在10欧姆附近,因此后匹配电路的设计并不难(原版作者给出的我复现不了,效果较差,因此我自己设计了一遍):
3.5、联合仿真
全部设计完了,联合仿真就行啦:
效果还可以,最上面已经放过结果了,这边多放一些图,比如饱和效率、增益、输出功率、回波损耗图等等:
4、文件解释(下载资源前必看)
4.1 环境
ADS版本:ADS2023
ADS依赖:VTB仿真平台(非必要),详见番外10:使用ADS对射频功率放大器进行非线性测试2(使用带宽20MHz的64QAM信号进行ACLR、EVM、CCDF测试)
依赖的ADS库文件:Cree公司的GAN库,CGH40010F管子、村田公司电容库(文件中包含,但可能要调整路径)
4.2 文件解释
AllCir文件夹下的All_Cir原理图:中途调试的PCB版图,无连接结构
AllCir文件夹下的All_Cir_v1原理图:最终的PCB版图
Input文件夹下的InputMatchCir原理图,是输入匹配的PCB版图EM仿真
Input文件夹下的InputMatchSub原理图,是输入匹配的原理图优化
Input文件夹下的InputMatchTest原理图,是输入匹配的版图结果的EM测试
Input文件夹下的InputPhaseShifter原理图,是输入匹配的相位延迟线的测试
Input文件夹下的SourcePull原理图,是源牵引的模板
InputRaw文件夹中是作者原来的输入匹配版图,我复现的效果不好就不额外赘述了。
Output文件夹下的FinalOutputMatch原理图,是后匹配电路的原理图仿真优化
Output文件夹下的OutputMatchCir原理图,是输出匹配的PCB版图EM仿真
Output文件夹下的OnputMatchSub原理图,是输出匹配的原理图的仿真(用的作者提供的参数)
Output文件夹下的OutputSubAuxiliary原理图,是对辅助功放后微带结构延迟的验证
Output文件夹下的OutputSubMain原理图,是对主功放后微带结构延迟的验证
Test文件夹中的DohertyAmp原理图,是DPA系统的版图联合仿真,但是所有部件没有放在一个板子上
Test文件夹中的DohertyAmpAll原理图,是DPA系统的版图联合仿真,所有部件放在一个板子上
Test文件夹中的DohertyAmpAllFinal原理图,是DPA系统的版图联合仿真,增加了实际焊接的焊盘等等
Test文件夹中的DohertyAmpAllFinal_v2原理图,是DPA系统的版图联合仿真,使用了实际的电阻电容模型
Test文件夹中的DohertyAmpRaw原理图,是作者提供原理图的复现,效果不好
Test文件夹中的DPD_Test原理图,使用VTB对连续信号进行测试,依赖于VTB仿真平台
Test文件夹中的HB1ToneGComp1swp_Doherty原理图,验证性能的原理图,非常重要!!!,选择验证哪一个版图联合仿真的性能,推荐验证最终的DohertyAmpAllFinal_v2的就行:
Test文件夹中的HB1TonePAE_Pswp_Doherty原理图,单频点的测试,一般不用。
Wilkinson文件夹是功分器的设计,比较简单就不说了
C_test是测试实际电容的谐振特性的,用于选择实际的电容
5、Matlab代码
% Single-Input Broadband Hybrid Doherty Power Amplifiers Design Relying
% on a Phase Sliding-Mode of the Load Modulation Scheme
clc
clear
close allf_arrary = 2.5:0.1:3.2;
K_va = [1.5, 1.75, 1.97, 2.15, 2.30, 2.40, 2.47, 2.51];
OBO_dB=8;
OBO=10^(OBO_dB/10);%% theta_b_cal
% EQN3
theta_b1= acosd(sqrt((OBO-K_va.^2)./(OBO+2*K_va+1)));
theta_b2= acosd(-sqrt((OBO-K_va.^2)./(OBO+2*K_va+1)));
theta_b3= -acosd(sqrt((OBO-K_va.^2)./(OBO+2*K_va+1)));
theta_b4= -acosd(-sqrt((OBO-K_va.^2)./(OBO+2*K_va+1)));
% cos计算,算出负值直接取反,cos(x)=cos(-x)
if theta_b1<0 theta_b1=-theta_b1;end
if theta_b2<0 theta_b2=-theta_b2;end
if theta_b3<0 theta_b3=-theta_b3;end
if theta_b4<0 theta_b4=-theta_b4;endfigure(3)
subplot(2,2,2);
plot(f_arrary,theta_b1,'-o','Color','r','LineWidth',3);hold on
plot(f_arrary,theta_b2,'-square','Color','magenta','LineWidth',3);hold on
plot(f_arrary,theta_b3,'-v','Color','b','LineWidth',1);hold on
plot(f_arrary,theta_b4,'-o','Color','black','LineWidth',1);
legend('\theta_{b1}(++)','\theta_{b2}(+-)','\theta_{b3}(-+)','\theta_{b4}(--)')
axis([min(f_arrary) max(f_arrary) 40 140])theta_p1=180-theta_b1;theta_p2=180-theta_b2;theta_p3=180-theta_b3;theta_p4=180-theta_b4;
subplot(2,2,1);
plot(f_arrary,theta_p1,'-o','Color','r','LineWidth',3);hold on
plot(f_arrary,theta_p2,'-square','Color','magenta','LineWidth',3);hold on
plot(f_arrary,theta_p3,'-v','Color','b','LineWidth',1);hold on
plot(f_arrary,theta_p4,'-o','Color','black','LineWidth',1);
legend('\theta_{p1}(++)','\theta_{p2}(+-)','\theta_{p3}(-+)','\theta_{p4}(--)')
axis([min(f_arrary) max(f_arrary) 40 140])%% theta_m_cal
% EQN4
theta_m1 = 180*atan((K_va+OBO)./(K_va+1).*tand(theta_b1))/pi;
theta_m2 = 180*atan((K_va+OBO)./(K_va+1).*tand(theta_b2))/pi;
theta_m3 = 180*atan((K_va+OBO)./(K_va+1).*tand(theta_b3))/pi;
theta_m4 = 180*atan((K_va+OBO)./(K_va+1).*tand(theta_b4))/pi;
% tan计算,算出负值直接加2pi
if theta_m1<0 theta_m1=theta_m1+360;end
if theta_m2<0 theta_m2=theta_m2+360;end
if theta_m3<0 theta_m3=theta_m3+360;end
if theta_m4<0 theta_m4=theta_m4+360;endsubplot(2,2,3);
plot(f_arrary,theta_m1,'-o','Color','r','LineWidth',2);hold on
plot(f_arrary,theta_m2,'-square','Color','magenta','LineWidth',2);hold on
plot(f_arrary,theta_m3,'-v','Color','b','LineWidth',1);hold on
plot(f_arrary,theta_m4,'-o','Color','black','LineWidth',1);
legend('\theta_{m1}(++)','\theta_{m2}(+-)','\theta_{m3}(-+)','\theta_{m4}(--)')
axis([min(f_arrary) max(f_arrary) 50 300])
%% theta_a_cal
% EQN8
% tan计算,算出来加pi结果不变
theta_a1 = atand(tand(theta_m1).*(K_va-OBO./K_va)./(K_va+OBO))+180;
theta_a2 = atand(tand(theta_m2).*(K_va-OBO./K_va)./(K_va+OBO))+180;
theta_a3 = atand(tand(theta_m3).*(K_va-OBO./K_va)./(K_va+OBO))+180;
theta_a4 = atand(tand(theta_m4).*(K_va-OBO./K_va)./(K_va+OBO))+180;subplot(2,2,4);
plot(f_arrary,theta_a1,'-o','Color','r','LineWidth',2);hold on
plot(f_arrary,theta_a2,'-square','Color','magenta','LineWidth',2);hold on
plot(f_arrary,theta_a3,'-v','Color','b','LineWidth',1);hold on
plot(f_arrary,theta_a4,'-o','Color','black','LineWidth',1);
legend('\theta_{a1}(++)','\theta_{a2}(+-)','\theta_{a3}(-+)','\theta_{a4}(--)')
axis([min(f_arrary) max(f_arrary) 120 240])%% Chapter2_2
Rmp=24.4;%Ropt
n_omega=K_va*(OBO-1)./(OBO+K_va);
RL=Rmp./(n_omega+1);
%使用并联的匹配关系计算Rap,RL=Rap//Rmp
Rap=1./(1./RL-1/Rmp);
%计算回退阻抗Rmb,先计算峰值开路微带线TL2阻抗
Zal_1=Rap.*(1./(1j*(tand(theta_a1))));
%TL2阻抗和RL并联
Zal_2=1./(1./Zal_1+1./RL);
%TL2和RL并联阻抗,与TL1阻抗变换到Rmb
Rmb=(Rmp*(Zal_2+Rmp*1j.*tand(theta_m1))./(Rmp+Zal_2.*1j.*tand(theta_m1)));
%Rmb虚部约等于0,直接忽略
Rmb=abs(Rmb);%EQN9_Cal
Z11=Rmb;
Z12=Rap*(OBO-1)./(K_va+1).*exp(-1j*theta_b1*2*pi/360);
Z21=Rap*(OBO-1)./(K_va+1).*exp(-1j*theta_b1*2*pi/360);
beta=(1+(OBO+K_va)./(K_va.*(K_va+1))+(K_va.^2-OBO)./(K_va.*(K_va+1)).*1j.*tand(theta_b1))./(1+1j.*tand(theta_b1));
Z22=Rap.*beta;
%计算输入信号的延迟,实际上就是两个微带线TL1、TL2的延迟
theta=abs(theta_a1-theta_m1);
%EQN10_Cal&Figure4 Draw;Im直接默认为1
Im=1;
Ia=1*exp(-1j*theta_b1*pi/180);
Vm=Z11*Im+Z12.*Ia;
Va=Z12*Im+Z22.*Ia;
Pout=0.5*real(Vm.*conj(Im))+0.5*real(Va.*conj(Ia));
delta_P1=10*log10(Pout/min(Pout));Rmp=24.4;%Ropt
% theta=180-theta_b1;
% OBO=(-K_va.*K_va-2*K_va.*cosd(theta).*cosd(theta)-cosd(theta).*cosd(theta))./(cosd(theta).*cosd(theta)-1)
n_omega=K_va*(OBO-1)./(OBO+K_va);
RL=Rmp./(n_omega+1);
%使用并联的匹配关系计算Rap,RL=Rap//Rmp
Rap=1./(1./RL-1/Rmp);
%计算回退阻抗Rmb,先计算峰值开路微带线TL2阻抗
Zal_1=Rap.*(1./(1j*(tand(theta_a4))));
%TL2阻抗和RL并联
Zal_2=1./(1./Zal_1+1./RL);
%TL2和RL并联阻抗,与TL1阻抗变换到Rmb
Rmb=(Rmp*(Zal_2+Rmp*1j.*tand(theta_m4))./(Rmp+Zal_2.*1j.*tand(theta_m4)));
%Rmb虚部约等于0,直接忽略
Rmb=abs(Rmb);Z11=Rmb;
Z12=Rap*(OBO-1)./(K_va+1).*exp(-1j*theta_p4*2*pi/360);
Z21=Rap*(OBO-1)./(K_va+1).*exp(-1j*theta_p4*2*pi/360);
beta=(1+(OBO+K_va)./(K_va.*(K_va+1))+(K_va.^2-OBO)./(K_va.*(K_va+1)).*1j.*tand(theta_p4))./(1+1j.*tand(theta_p4));
Z22=Rap.*beta;
Ia=1*exp(-1j*theta_p4*pi/180);
Vm=Z11*Im+Z12.*Ia;
Va=Z12*Im+Z22.*Ia;
Pout=0.5*real(Vm.*conj(Im))+0.5*real(Va.*conj(Ia));
delta_P2=10*log10(Pout/min(Pout));figure(4)
plot(K_va,delta_P1,'-o','Color','magenta','LineWidth',2);hold on
axis([min(K_va) max(K_va) 0 3.5])
plot(K_va,delta_P2,'-v','Color','black','LineWidth',2);hold on
legend('\theta_{b1}(++)','\theta_{p4}(--)')