【负荷预测】基于蚂蚁优化算法的BP神经网络在负荷预测中的应用研究（Matlab完整代码实现）

news2024/11/15 9:05:14

目录

0 知识回顾

1 ACO-BP算法

2 ACO-BP算法基本思路

3 具体步骤

4 Matlab代码实现

5 运行结果

6 参考文献 

7 写在最后 

0 知识回顾

智能优化算法—蚁群算法（Python实现）

1 ACO-BP算法


传统的BP神经网络训练采用的是误差反向传播学习算法,它的优化目标函数相对复杂,较容易出现陷人局部最优、收敛速度慢等问题[6]。由于BP神经网络的训练算法实质上是对其网络权值和阈值进行迭代调整,因此用蚁群优化算法替代BP算法完成对神经网络权阈值的迭代调整,并最终完成神经网络的训练。

2 ACO-BP算法基本思路


蚁群算法解决优化问题的基本思路为:用蚂蚁的行走路径表示待优化问题的可行解,整个蚂蚁群体的所有路径构成待优化问题的解空间。路径较短的蚂蚁释放较多的信息素。经过一定时间,信息素浓度在较短的路径上累计较高,所有选择此路径的蚂蚁也逐渐增多,最终,整个蚁群会在正反馈的作用下集中在最佳路径上,此时对应的便是待优化问题的最优解。
首先根据权值和阈值的取值池间,将的E义以刘你以S个等长区间,即每个区间的长度被作s等分,将区间的临界值或选择区间中的随机值作为候选值。确定参数个数n,包括网络中所有的权阈值。每个参数Pi(i=1,2,…,n)对应个有S个元素的集合l, ,这些元素为Pi的可能取值。

3 具体步骤


1）参数初始化:将所有权值和阈值进行S等分,所有区间初始信息素О,信息素残留系教入,1后尽系代i,前区间信息表为Tabu,最大迭代次数C,网络全局误差E，最大学习次数N;
2)权值和阈值选择:蚁群m只蚂蚁,对于蚂蚁k依据概率公式(2)的寻路规则进行选择节点所在区间,蚁群迭代一次则完成一次解的构造:

            

式中\tau _{j}(I_{pi})为集合I_{pi}中第j个元素的信息素值;
3）蚁群寻优判断:蚁群迭代一次得到的构造解,则是当前迭代后得到误差最小的一组解,计算误差Ec ,判断是否达到蚁群要求,若是则转到4),否则转到5);
4)网络训练:将蚁群迭代得到的最优构造解,作为初始权值和阈值,选取数据集对网络进行训练,直到满足结束条件即最大学习次数,完成学习。否则,继续学习;
5）更新信息素:根据式(3)、式( 4)、式(5)对所有区间信息素全局更新,并重置信息表:

              

6)蚁群遍历:重复步骤2)到步骤3)。其算法如图1所示。
           

                                图1 ACO-BP算法流程图

4 Matlab代码实现

主函数

%% 清空环境变量
clc
clear
close all
format compact 
%% 网络结构建立
%% 清空环境变量
clc
clear
close all
format compact 
%% 网络结构建立
%读取数据
data=xlsread('天气_电量_数据.xlsx','C12:J70');%前7列为每个时刻的发电量 最后列为天气

for i=1:58
    input(i,:)=[data(i,:) data(i+1,end)];
    output(i,:)=data(i+1,1:7);
end

%% 节点个数
inputnum=9;%输入 前一天7个时刻的电量+前一天的天气+预测日的天气
hiddennum=5;
outputnum=7;%预测日7个时刻的发电量

%% 训练数据和预测数据 最后一天用来测试  前面的都拿来训练
input_train=input(1:57,:)';
input_test=input(58,:)';
output_train=output(1:57,:)';
output_test=output(58,:)';

%选连样本输入输出数据归一化
[inputn,inputps]=mapminmax(input_train);
[outputn,outputps]=mapminmax(output_train);
inputn_test=mapminmax('apply',input_test,inputps);
%%
%构建网络
net=newff(inputn,outputn,hiddennum);
%寻优
[bestnest,trace]=antforelm(inputnum,hiddennum,outputnum,net,inputn,outputn);
figure
plot(trace)
title('适应度曲线')
xlabel('迭代数')
ylabel('适应度值')


%% 把最优初始阀值权值赋予网络预测
x=bestnest;
% 用CS优化的BP网络进行值预测
w1=x(1:inputnum*hiddennum);
B1=x(inputnum*hiddennum+1:inputnum*hiddennum+hiddennum);
w2=x(inputnum*hiddennum+hiddennum+1:inputnum*hiddennum+hiddennum+hiddennum*outputnum);
B2=x(inputnum*hiddennum+hiddennum+hiddennum*outputnum+1:inputnum*hiddennum+hiddennum+hiddennum*outputnum+outputnum);

net.iw{1,1}=reshape(w1,hiddennum,inputnum);
net.lw{2,1}=reshape(w2,outputnum,hiddennum);
net.b{1}=reshape(B1,hiddennum,1);
net.b{2}=B2';

%% BP网络训练
%网络进化参数
net.trainParam.epochs=200;
net.trainParam.lr=0.1;
%net.trainParam.goal=0.00001;

%网络训练
[net,per2]=train(net,inputn,outputn);

%% BP网络预测
%数据归一化
inputn_test=mapminmax('apply',input_test,inputps);
an=sim(net,inputn_test);
test_simu=mapminmax('reverse',an,outputps);
error=test_simu-output_test;
%%
figure
a1=output_test;
a2=test_simu;
plot(a1,'*-');hold on
plot(a2,'O-')
title('各时刻发电量实际值与预测值')
xlabel('')
legend('原始数据','bp预测数据')
set(gca,'XTick',1:7,...                                    
        'XTickLabel',{'9:00','10:00','11:00','12:00','13:00','14:00','15:00'},...
        'TickLength',[0 0]);
grid on
ylabel('发电量（KW）')

%% You can replace the following by your own functions

function error = fun(x,inputnum,hiddennum,outputnum,net,inputn,outputn)
%该函数用来计算适应度值
%x          input     个体
%inputnum   input     输入层节点数
%outputnum  input     隐含层节点数
%net        input     网络
%inputn     input     训练输入数据
%outputn    input     训练输出数据
%error      output    个体适应度值

%提取
w1=x(1:inputnum*hiddennum);
B1=x(inputnum*hiddennum+1:inputnum*hiddennum+hiddennum);
w2=x(inputnum*hiddennum+hiddennum+1:inputnum*hiddennum+hiddennum+hiddennum*outputnum);
B2=x(inputnum*hiddennum+hiddennum+hiddennum*outputnum+1:inputnum*hiddennum+hiddennum+hiddennum*outputnum+outputnum);


%网络进化参数
net.trainParam.epochs=20;
net.trainParam.lr=0.1;
net.trainParam.goal=0.00001;
net.trainParam.show=100;
net.trainParam.showWindow=0;
 
%网络权值赋值
net.iw{1,1}=reshape(w1,hiddennum,inputnum);
net.lw{2,1}=reshape(w2,outputnum,hiddennum);
net.b{1}=reshape(B1,hiddennum,1);
net.b{2}=B2';

% %网络训练
% net=train(net,inputn,outputn);
an=sim(net,inputn);
% an-outputn
error0=sum(abs(an-outputn));
error=sum(error0);%以均方差作为适应度函数

%% Replace some nests by constructing new solutions/nests
function new_nest=empty_nests(nest,Lb,Ub,pa)
% A fraction of worse nests are discovered with a probability pa
n=size(nest,1);
% Discovered or not -- a status vector
K=rand(size(nest))>pa;

% In the real world, if a cuckoo's egg is very similar to a host's eggs, then 
% this cuckoo's egg is less likely to be discovered, thus the fitness should 
% be related to the difference in solutions.  Therefore, it is a good idea 
% to do a random walk in a biased way with some random step sizes.  
%% New solution by biased/selective random walks
stepsize=rand*(nest(randperm(n),:)-nest(randperm(n),:));
new_nest=nest+stepsize.*K;
for j=1:size(new_nest,1)
    s=new_nest(j,:);
  new_nest(j,:)=simplebounds(s,Lb,Ub);  
end

function [y,trace]=antforelm(inputnum,hiddennum,outputnum,net,inputn_train,label_train);%蚁群算法%%%%%%%%%%%%%%%%%%%%蚁群算法求函数极值%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%初始化%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
m=5;                    %蚂蚁个数
G_max=100;               %最大迭代次数
Rho=0.5;                 %信息素蒸发系数
P0=0.5;                  %转移概率常数
XMAX= 1;                 %搜索变量x最大值
XMIN=-1;                %搜索变量x最小值
d=inputnum*hiddennum+hiddennum+hiddennum*outputnum+outputnum;    
%%%%%%%%%%%%%%%%%随机设置蚂蚁初始位置%%%%%%%%%%%%%%%%%%%%%%
for i=1:m
    X(i,:)=(XMIN+(XMAX-XMIN).*rand(1,d));
    Tau(i)=fun(X(i,:),inputnum,hiddennum,outputnum,net,inputn_train,label_train); 
end

bestfitness=inf;
bestfitness_position=inf*ones(1,d);

step=0.1;                %局部搜索步长
for NC=1:G_max
    NC
    lamda=1/NC;
    [Tau_best,BestIndex]=min(Tau);
    %%%%%%%%%%%%%%%%%%计算状态转移概率%%%%%%%%%%%%%%%%%%%%
    for i=1:m
        P(NC,i)=(Tau(BestIndex)-Tau(i))/Tau(BestIndex);
    end
    %%%%%%%%%%%%%%%%%%%%%%位置更新%%%%%%%%%%%%%%%%%%%%%%%%
    for i=1:m
        fun_i=fun(X(i,:),inputnum,hiddennum,outputnum,net,inputn_train,label_train);
           %%%%%%%%%%%%%%%%%局部搜索%%%%%%%%%%%%%%%%%%%%%%
        if P(NC,i)<P0
            temp1=X(i,:)+(rand(1,d))*step*lamda;
           
        else
            %%%%%%%%%%%%%%%%全局搜索%%%%%%%%%%%%%%%%%%%%%%%
             temp1=X(i,:)+(XMAX-XMIN)*(rand(1,d));
        end
        %%%%%%%%%%%%%%%%%%%%%边界处理%%%%%%%%%%%%%%%%%%%%%%%
        for j=1:d
            if temp1(j)<XMIN
                temp1(j)=rand;
            end
            if temp1(j)>XMAX
                temp1(j)=rand;
            end
        end
        fun_temp=fun(temp1,inputnum,hiddennum,outputnum,net,inputn_train,label_train);
        %%%%%%%%%%%%%%%%%%蚂蚁判断是否移动%%%%%%%%%%%%%%%%%%
        if fun_temp<fun_i
            X(i,:)=temp1;
            Tau(i)=fun_temp;
        end
    end
    %%%%%%%%%%%%%%%%%%%%%%%更新信息素%%%%%%%%%%%%%%%%%%%%%%%
    for i=1:m
        Tau(i)=(1-Rho)*Tau(i)+Tau(i); 
    end
    [value,index]=min(Tau);
    %%
    if value<bestfitness
        bestfitness=value;
        bestfitness_position=X(index,:);
    end
    trace(NC)=bestfitness;  
    
end
y=bestfitness_position;                           %最优变量

5 运行结果

 

 

6 参考文献 

[1]陈智雨,陆金桂.基于ACO-BP神经网络的光伏系统发电功率预测[J].机械制造与自动化,2020,49(01):173-175+187.DOI:10.19344/j.cnki.issn1671-5276.2020.01.047.

7 写在最后 

部分理论引用网络文献，如有侵权请联系删除。

%% 清空环境变量
clc
clear
close all
format compact 
%% 网络结构建立
%% 清空环境变量
clc
clear
close all
format compact 
%% 网络结构建立
%读取数据
data=xlsread('天气_电量_数据.xlsx','C12:J70');%前7列为每个时刻的发电量 最后列为天气

for i=1:58
    input(i,:)=[data(i,:) data(i+1,end)];
    output(i,:)=data(i+1,1:7);
end

%% 节点个数
inputnum=9;%输入 前一天7个时刻的电量+前一天的天气+预测日的天气
hiddennum=5;
outputnum=7;%预测日7个时刻的发电量

%% 训练数据和预测数据 最后一天用来测试  前面的都拿来训练
input_train=input(1:57,:)';
input_test=input(58,:)';
output_train=output(1:57,:)';
output_test=output(58,:)';

%选连样本输入输出数据归一化
[inputn,inputps]=mapminmax(input_train);
[outputn,outputps]=mapminmax(output_train);
inputn_test=mapminmax('apply',input_test,inputps);
%%
%构建网络
net=newff(inputn,outputn,hiddennum);
%寻优
[bestnest,trace]=antforelm(inputnum,hiddennum,outputnum,net,inputn,outputn);
figure
plot(trace)
title('适应度曲线')
xlabel('迭代数')
ylabel('适应度值')


%% 把最优初始阀值权值赋予网络预测
x=bestnest;
% 用CS优化的BP网络进行值预测
w1=x(1:inputnum*hiddennum);
B1=x(inputnum*hiddennum+1:inputnum*hiddennum+hiddennum);
w2=x(inputnum*hiddennum+hiddennum+1:inputnum*hiddennum+hiddennum+hiddennum*outputnum);
B2=x(inputnum*hiddennum+hiddennum+hiddennum*outputnum+1:inputnum*hiddennum+hiddennum+hiddennum*outputnum+outputnum);

net.iw{1,1}=reshape(w1,hiddennum,inputnum);
net.lw{2,1}=reshape(w2,outputnum,hiddennum);
net.b{1}=reshape(B1,hiddennum,1);
net.b{2}=B2';

%% BP网络训练
%网络进化参数
net.trainParam.epochs=200;
net.trainParam.lr=0.1;
%net.trainParam.goal=0.00001;

%网络训练
[net,per2]=train(net,inputn,outputn);

%% BP网络预测
%数据归一化
inputn_test=mapminmax('apply',input_test,inputps);
an=sim(net,inputn_test);
test_simu=mapminmax('reverse',an,outputps);
error=test_simu-output_test;
%%
figure
a1=output_test;
a2=test_simu;
plot(a1,'*-');hold on
plot(a2,'O-')
title('各时刻发电量实际值与预测值')
xlabel('')
legend('原始数据','bp预测数据')
set(gca,'XTick',1:7,...                                    
        'XTickLabel',{'9:00','10:00','11:00','12:00','13:00','14:00','15:00'},...
        'TickLength',[0 0]);
grid on
ylabel('发电量（KW）')
 

本文来自互联网用户投稿，该文观点仅代表作者本人，不代表本站立场。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如若转载，请注明出处：http://www.coloradmin.cn/o/56235.html

如若内容造成侵权/违法违规/事实不符，请联系多彩编程网进行投诉反馈，一经查实，立即删除！

相关文章

2023最新SSM计算机毕业设计选题大全（附源码+LW）之java高校教室申请管理系统yf80k

2023最新SSM计算机毕业设计选题大全（附源码+LW）之java高校教室申请管理系统yf80k

计算机毕业设计的编程真的不会, 应该怎么办啊, 平时学了3年都没学懂&#xff0c;然而还有一个月就要答辩了&#xff0c;一点东西都没做出来&#xff0c;不知道该怎么办了&#xff1f;好迷茫&#xff0c;如果毕不到业怎么办, 有没有快速的学习方法 毕设毕设&#xff0c;最终的是…
阅读更多...
【C++面向对象程序设计】CH7 输入输出流

【C++面向对象程序设计】CH7 输入输出流

目录 一、前言 二、标准输出流 1.流对象 &#xff08;1&#xff09;cout流对象 &#xff08;2&#xff09;cerr流对象 &#xff08;3&#xff09;【例7.1】编写程序&#xff0c;从键盘输入a&#xff0c;b&#xff0c;c的值求解一元二次方程。如果a0或判别式的值<0&…
阅读更多...
拾壹博客拆解改造，页面元素替换（二）

拾壹博客拆解改造，页面元素替换（二）

页面元素替换 首先要做的当然是换成自己风格的站名和内容啦。 1、网站配置 跟踪前端代码后发现配置是来自后端接口&#xff0c;想着既然入库了&#xff0c;那应该有对应的管理页面吧&#xff0c;果然找到了&#xff0c;就是…演示账号不允许操作&#xff01;那么接下来要干的…
阅读更多...
stm32f767的fft

stm32f767的fft

仅作自己笔记用 1&#xff0c;FFT函数调用基础知识 采样得到的数字信号&#xff0c;就可以做FFT变换了。N个采样点&#xff0c;经过FFT之后&#xff0c;就可以得到N个点的FFT结果。为了方便进行FFT运算&#xff0c;通常N取2的整数次方。 假设采样频率为Fs&#xff0c;信号频率…
阅读更多...
基于Java环境下的高校跳蚤市场商城系统

基于Java环境下的高校跳蚤市场商城系统

目 录 摘 要 I Abstract II 1绪论 1 1.1 课题背景 1 1.2 目的和意义 1 1.3 研究现状 2 1.4 研究主要内容 3 2开发平台与技术的介绍 4 2.1 Eclipse简介 4 2.2 Java EE简介 4 2.2.1 Java EE概念 4 2.2.2 Java EE运行模式 4 2.3 Jsp技术简介 5 2.4 Struts 2框架简介 5 2.5 MySQL简…
阅读更多...
[附源码]计算机毕业设计springboot新冠疫苗接种预约系统

[附源码]计算机毕业设计springboot新冠疫苗接种预约系统

项目运行 环境配置&#xff1a; Jdk1.8 Tomcat7.0 Mysql HBuilderX&#xff08;Webstorm也行&#xff09; Eclispe&#xff08;IntelliJ IDEA,Eclispe,MyEclispe,Sts都支持&#xff09;。 项目技术&#xff1a; SSM mybatis Maven Vue 等等组成&#xff0c;B/S模式 M…
阅读更多...
Docker最新超详细教程——Docker创建运行Redis并挂载

Docker最新超详细教程——Docker创建运行Redis并挂载

Docker最新超详细教程——Docker创建运行Redis并挂载 Docker官网关于Redis的描述 redis - Official Image | Docker Hubhttps://hub.docker.com/_/redis 一、拉取Redis镜像 docker pull <镜像名称>:<版本号> docker pull redis:6.2.27 首先我们要在Docker上获取…
阅读更多...
2023最新SSM计算机毕业设计选题大全（附源码+LW）之java高校教学过程管理系统34085

2023最新SSM计算机毕业设计选题大全（附源码+LW）之java高校教学过程管理系统34085

现在毕设刚开始。时间还有很多&#xff0c;可以从头开始学也可以。毕设其实不难&#xff0c;难的是我们懒散到这种时候再去静下心学。能自己独立完成尽量自己独立完成。相信你看过很多上面回答的&#xff0c;都不建议去某宝。毕竟这一行参差不齐哈。能找到靠谱的也不容易。近期…
阅读更多...
Unity Debug的简单封装

Unity Debug的简单封装

对Unity Debug的简单封装 使用前提&#xff1a; Project Settings-Player-Other Settings-Script Define Symbols添加 EnableLog&#xff0c;点击Apply 测试代码&#xff1a; using MTools.Debuger; using UnityEngine;public class NewBehaviourScript : MonoBehaviour {p…
阅读更多...
基于PHP+MySQL家庭医生签约预约诊疗管理信息系统

基于PHP+MySQL家庭医生签约预约诊疗管理信息系统

随着时代的发展,人们对医疗方面的要求也越来越高,也是人们更希望通过家庭医生来对自己提供所需的医疗服务,从而享受更加个性化的医疗服务,为此我们开发了本家庭医生签约预约诊疗管理信息系统,通过本系统患者可以享有签约,预约,和诊疗等一系类的服务。 本系统是一个家庭医生签约…
阅读更多...
[附源码]Python计算机毕业设计Django基于微信小程序的网络办公系统

[附源码]Python计算机毕业设计Django基于微信小程序的网络办公系统

项目运行 环境配置&#xff1a; Pychram社区版 python3.7.7 Mysql5.7 HBuilderXlist pipNavicat11Djangonodejs。 项目技术&#xff1a; django python Vue 等等组成&#xff0c;B/S模式 pychram管理等等。 环境需要 1.运行环境&#xff1a;最好是python3.7.7&#xff0c;…
阅读更多...
qmake 参数

qmake 参数

E:\workspace\QtWork\qmake\option.cpp:Option::init() -project 设置qmake生成工程相关文件&#xff0c;如果用qt creator开发的话这个命令参数基本用不到。 -prl 设置qmake生成prl文件。 -set 设置自定义属性&#xff0c;会存放到注册表中。具体参考属性 -unset 取消自定义…
阅读更多...
网页JS自动化脚本(六)在特定位置添加元素

网页JS自动化脚本(六)在特定位置添加元素

在某元素后插入元素 我们这一次在按钮元素后面复制一个一模一样的按钮,所以分为几个步骤,先新建一个一样的元素,然后把相同中的属性赋值给它,再插入到合适的位置,最后再稍微修改一下外观样式即可 首先新建一个input元素,看一下效果 window.onloadfunction(){var theElementdo…
阅读更多...
[附源码]Python计算机毕业设计Django基于人脸识别的社区防疫管理系统

[附源码]Python计算机毕业设计Django基于人脸识别的社区防疫管理系统

项目运行 环境配置&#xff1a; Pychram社区版 python3.7.7 Mysql5.7 HBuilderXlist pipNavicat11Djangonodejs。 项目技术&#xff1a; django python Vue 等等组成&#xff0c;B/S模式 pychram管理等等。 环境需要 1.运行环境&#xff1a;最好是python3.7.7&#xff0c;…
阅读更多...
Pytest自动化测试框架---（单元测试框架）

Pytest自动化测试框架---（单元测试框架）

unittest是python自带的单元测试框架&#xff0c;它封装好了一些校验返回的结果方法和一些用例执行前的初始化操作&#xff0c;使得单元测试易于开展&#xff0c;因为它的易用性&#xff0c;很多同学也拿它来做功能测试和接口测试&#xff0c;只需简单开发一些功能&#xff08;…
阅读更多...
阿里云新用户活动：云服务器ECS 新购、升级报价出炉了！

阿里云新用户活动：云服务器ECS 新购、升级报价出炉了！

阿里云新人特惠&#xff0c;阿里云新用户新购升级立享满减&#xff0c;新购升级云服务器ECS &#xff0c;购买热门产品 s6/u1/c6/g6/r6/c7/g7/r7指定配置&#xff0c;可享折上折&#xff01;从未购买过云服务器ECS或者轻量应用服务器的用户一次性可领取3张优惠券。优惠券适用于…
阅读更多...
深入浅出Seata的AT模式

深入浅出Seata的AT模式

单个掉队&#xff0c;导致集体被动摆烂&#xff1b; 一、业务背景 在分布式架构中&#xff0c;事务管理是个无法避开的复杂问题&#xff0c;虽然有多种解决方案&#xff0c;但是需要根据业务去选择合适的&#xff1b; 从个人最近几年的实践经验来看&#xff0c;Seata组件的AT…
阅读更多...
【JUC】SpringBoot使用线程池的两种方式 注解和注入

【JUC】SpringBoot使用线程池的两种方式 注解和注入

学习笔记一、ThreadPoolTaskExecutor与ThreadPoolExecutor的区别二、编写配置文件ThreadPoolConfig二、编写Controller三、编写Service3.1、注解3.1、注入一、ThreadPoolTaskExecutor与ThreadPoolExecutor的区别 ThreadPoolExecutor 是JDK自1.5添加的线程池。 ThreadPoolTaskE…
阅读更多...
排序算法:插入排序,选择排序,冒泡排序

排序算法:插入排序,选择排序,冒泡排序

插入排序 一般来说&#xff0c;插入排序都采用in-place在数组上实现。具体算法描述如下&#xff1a; 步骤1: 从第一个元素开始&#xff0c;该元素可以认为已经被排序&#xff1b; 步骤2: 取出下一个元素&#xff0c;在已经排序的元素序列中从后向前扫描&#xff1b; 步骤3: 如…
阅读更多...
[附源码]计算机毕业设计-中国传统手工艺销售平台Springboot程序

[附源码]计算机毕业设计-中国传统手工艺销售平台Springboot程序

项目运行 环境配置&#xff1a; Jdk1.8 Tomcat7.0 Mysql HBuilderX&#xff08;Webstorm也行&#xff09; Eclispe&#xff08;IntelliJ IDEA,Eclispe,MyEclispe,Sts都支持&#xff09;。 项目技术&#xff1a; SSM mybatis Maven Vue 等等组成&#xff0c;B/S模式 M…
阅读更多...
推荐文章
最新文章

Copyright @ 2022~2023