基于Qlearning强化学习的倒立摆控制系统matlab仿真1.算法描述2.仿真效果预览3.MATLAB部分代码预览4.完整MATLAB程序

131 阅读 0 评论 87 点赞

我是靠谱客的博主谨慎小蘑菇，这篇文章主要介绍基于Qlearning强化学习的倒立摆控制系统matlab仿真1.算法描述2.仿真效果预览3.MATLAB部分代码预览4.完整MATLAB程序，现在分享给大家，希望可以做个参考。

1.算法描述

2.仿真效果预览

3.MATLAB部分代码预览

4.完整MATLAB程序

1.算法描述

强化学习通常包括两个实体agent和environment。两个实体的交互如下，在environment的statestst下，agent采取actionatat进而得到rewardrtrt 并进入statest+1st+1。Q-learning的核心是Q-table。Q-table的行和列分别表示state和action的值，Q-table的值Q(s,a)Q(s,a)衡量当前states采取actiona到底有多好。

在每一时刻，智能体观测环境的当下状态并选择一个动作，这会导致环境转移到一个新的状态，与此同时环境会返回给智能体一个奖励，该奖励反映了动作所导致的结果。在倒立摆任务中，每一个时间步的奖励均为+1，但是一旦小车偏离中心超过4.8个单位或者杆的倾斜超过15度，任务就会终止。因此，我们的目标是使得该任务能够尽可能地运行得更久，以便获得更多的收益。原始倒立摆任务中，智能体的输入包含4个实数（位置，速度等），但实际上，神经网络可以直接通过观察场景来完成任务，所以我们可以直接使用以小车为中心的屏幕补丁作为输入。严格来说，我们设计的状态是当前屏幕补丁与上一个屏幕补丁的差值，这使得智能体能够从一张图像中推断出杆的速度。

为了训练DQN，我们将使用经验回放池（experience replay memory）来存储智能体所观测到的环境状态转移情况，在之后的训练中我们可以充分利用这些数据。通过对经验回放池中的数据进行随机采样，组成一个批次的转移情况是互不相关（decorrelated）的，这极大地提升了DQN训练的性能和稳定性。

主要步骤如下：

采样得到一个批次的样本，将这些样本对应的张量连接成一个单独的张量；
分别利用策略Q网络与目标Q网络计算与Q(st,at)与V(st+1)=maxaQ(st+1,a) ，利用它们计算损失函数.。另外，如果 s 为终止状态，则令 V(s)=0
更新Q网络参数。目标Q网络的参数每隔一段时间从主Q网络处固定而来，在本例中，我们在每个episode更新一次目标Q网络。

2.仿真效果预览

matlab2022a仿真结果如下：

3.MATLAB部分代码预览

复制代码

for trial=1:MaxTr,      %外部循环开始 
    count=0; 
    failure=0; 
    failReason=0; 
    lfts = 1; 
    newSt = inistate; 
    inputs = newSt./NF; 
    lc = Initlc; 
    la = Initla; 
         
    xhist=newSt; 
    
    %计算newAction 
    ha = inputs*wa1; 
    g = (1 - exp(-ha))./(1 + exp(-ha)); 
    va = g*wa2; 
    newAction = (1 - exp(-va))./(1 + exp(-va)); 
    %计算J 
    inp=[inputs newAction]; 
    qc=inp*wc1; 
    p = (1 - exp(-qc))./(1 + exp(-qc)); 
    J=p*wc2; 
     
    Jprev = J; 
    
    while(lfts<Tit),        %内部循环开始 
             
         if (rem(lfts,500)==0), 
            disp(['It is ' int2str(lfts) ' time steps now......']); 
         end 
          
         %生成控制信号 
                if (newAction >= 0) 
                    sgnf = 1; 
                else 
                    sgnf = -1; 
                end 
                u = Mag*sgnf;		%bang-bang control 
 
		    %Plug in the model 
            [T,Xf]=ode45('cartpole_model',[0 tstep],newSt,[],u); 
	        a=size(Xf); 
	        newSt=Xf(a(1),:);         
            inputs=newSt./NF;	%input normalization    
             
                %计算newAction 
                ha = inputs*wa1; 
                g = (1 - exp(-ha))./(1 + exp(-ha)); 
                va = g*wa2; 
                newAction = (1 - exp(-va))./(1 + exp(-va)); 
                %calculate new J     
                inp=[inputs newAction]; 
                qc=inp*wc1; 
                p = (1 - exp(-qc))./(1 + exp(-qc)); 
                J=p*wc2; 
                 
                xhist=[xhist;newSt]; 
             
		    %%===========================================================%% 
		    %%求取强化信号r(t),即reinf                                   %% 
		    %%===========================================================%% 
          
            if (abs(newSt(1)) > FailTheta) 
                reinf = 1; 
                failure = 1; 
                failReason = 1; 
            elseif (abs(newSt(3)) > Boundary) 
                reinf = 1; 
                failure = 1; 
                failReason = 2; 
            else 
                reinf = 0; 
            end 
     
    	    %%================================%% 
    	    %% learning rate update scheme    %% 
    	    %%================================%% 
            
            if (rem(lfts,5)==0) 
                lc = lc - 0.05; 
                la = la - 0.05; 
            end 
           
      	    if (lc<0.01) 
                lc=0.005; 
            end 
           
            if (la<0.01) 
                la=0.005; 
            end 
           
            %%================================================%% 
            %% internal weights updating cycles for critnet   %% 
            %%================================================%% 
                    
            cyc = 0; 
            ecrit = alpha*J-(Jprev-reinf); 
		    Ec = 0.5 * ecrit^2; 
            while (Ec>Tc & cyc<=Ncrit), 
                    gradEcJ=alpha*ecrit; 
    				%----for the first layer(input to hidden layer)----------- 
                    gradqwc1 = [inputs'; newAction]; 
                    for i=1:N_Hidden, 
                        gradJp = wc2(i); 
                        gradpq = 0.5*(1-p(i)^2); 
                        wc1(:,i) = wc1(:,i) - lc*gradEcJ*gradJp*gradpq*gradqwc1; 
                    end        
                    %----for the second layer(hidden layer to output)----------- 
                    gradJwc2=p'; 
                    wc2 = wc2- lc*gradEcJ*gradJwc2; 
                    %----compute new  J---- 
                    inp=[inputs newAction]; 
                    qc=inp*wc1; 
                    p = (1 - exp(-qc))./(1 + exp(-qc)); 
                    J=p*wc2; 
 
                cyc = cyc +1; 
                ecrit = alpha*J-(Jprev-reinf); 
                Ec = 0.5 * ecrit^2; 
            end                                 % end of "while (Ec>0.05 & cyc<=Ncrit)" 
             
            %normalization weights for critical network 
                if (max(max(abs(wc1)))>1.5) 
                    wc1=wc1/max(max(abs(wc1))); 
                end 
                if max(max(abs(wc2)))>1.5 
                    wc2=wc2/max(max(abs(wc2))); 
                end 
                      
            %%=============================================%% 
            %% internal weights updating cycles for actnet %% 
            %%=============================================%%                 
            cyc = 0;             
            eact = J - Uc; 
            Ea = 0.5*eact^2; 
            while (Ea>Ta & cyc<=Nact), 
                graduv = 0.5*(1-newAction^2);              
                gradEaJ = eact; 
                gradJu = 0; 
                   for i=1:N_Hidden, 
                       gradJu = gradJu + wc2(i)*0.5*(1-p(i)^2)*wc1(WC_Inputs,i); 
                   end 
                    %----for the first layer(input to hidden layer)----------- 
                    for (i=1:N_Hidden), 
                        gradvg = wa2(i); 
                        gradgh = 0.5*(1-g(i)^2); 
                        gradhwa1 = inputs'; 
                        wa1(:,i)=wa1(:,i)-la*gradEaJ*gradJu*graduv*gradvg*gradgh*gradhwa1; 
                    end 
                    %----for the second layer(hidden layer to output)----------- 
                    gradvwa2 = g'; 
                    wa2=wa2-la*gradEaJ*gradJu*graduv*gradvwa2; 
                    %----compute new J and newAction-------  
                    ha = inputs*wa1; 
                    g = (1 - exp(-ha))./(1 + exp(-ha)); 
                    va = g*wa2; 
                    newAction = (1 - exp(-va))./(1 + exp(-va)); 
                     
                    inp=[inputs newAction]; 
                    qc=inp*wc1; 
                    p = (1 - exp(-qc))./(1 + exp(-qc)); 
                    J=p*wc2; 
                         
                cyc = cyc+1; 
                eact = J - Uc; 
                Ea = 0.5*eact^2;   
            end                       %end of "while (Ea>Ta & cyc<=Nact)" 
            
            if ~failure 
                Jprev=J; 
            else 
                break;                %another trial 即跳出“while(lfts<Tit),” 
            end 
            lfts=lfts+1; 
        end                           %end of "while(lfts<Tit)" 结束内部循环 
           
        msgstr1=['Trial # ' int2str(trial) ' has  ' int2str(lfts) ' time steps.']; 
   	    msgstr21=['Trial # ' int2str(trial) ' has successfully balanced for at least ']; 
   	    msgstr22=[msgstr21 int2str(lfts) ' time steps ']; 
A_027

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
 
for trial=1:MaxTr,      %外部循环开始 
    count=0; 
    failure=0; 
    failReason=0; 
    lfts = 1; 
    newSt = inistate; 
    inputs = newSt./NF; 
    lc = Initlc; 
    la = Initla; 
         
    xhist=newSt; 
    
    %计算newAction 
    ha = inputs*wa1; 
    g = (1 - exp(-ha))./(1 + exp(-ha)); 
    va = g*wa2; 
    newAction = (1 - exp(-va))./(1 + exp(-va)); 
    %计算J 
    inp=[inputs newAction]; 
    qc=inp*wc1; 
    p = (1 - exp(-qc))./(1 + exp(-qc)); 
    J=p*wc2; 
     
    Jprev = J; 
    
    while(lfts<Tit),        %内部循环开始 
             
         if (rem(lfts,500)==0), 
            disp(['It is ' int2str(lfts) ' time steps now......']); 
         end 
          
         %生成控制信号 
                if (newAction >= 0) 
                    sgnf = 1; 
                else 
                    sgnf = -1; 
                end 
                u = Mag*sgnf;		%bang-bang control 
 
		    %Plug in the model 
            [T,Xf]=ode45('cartpole_model',[0 tstep],newSt,[],u); 
	        a=size(Xf); 
	        newSt=Xf(a(1),:);         
            inputs=newSt./NF;	%input normalization    
             
                %计算newAction 
                ha = inputs*wa1; 
                g = (1 - exp(-ha))./(1 + exp(-ha)); 
                va = g*wa2; 
                newAction = (1 - exp(-va))./(1 + exp(-va)); 
                %calculate new J     
                inp=[inputs newAction]; 
                qc=inp*wc1; 
                p = (1 - exp(-qc))./(1 + exp(-qc)); 
                J=p*wc2; 
                 
                xhist=[xhist;newSt]; 
             
		    %%===========================================================%% 
		    %%求取强化信号r(t),即reinf                                   %% 
		    %%===========================================================%% 
          
            if (abs(newSt(1)) > FailTheta) 
                reinf = 1; 
                failure = 1; 
                failReason = 1; 
            elseif (abs(newSt(3)) > Boundary) 
                reinf = 1; 
                failure = 1; 
                failReason = 2; 
            else 
                reinf = 0; 
            end 
     
    	    %%================================%% 
    	    %% learning rate update scheme    %% 
    	    %%================================%% 
            
            if (rem(lfts,5)==0) 
                lc = lc - 0.05; 
                la = la - 0.05; 
            end 
           
      	    if (lc<0.01) 
                lc=0.005; 
            end 
           
            if (la<0.01) 
                la=0.005; 
            end 
           
            %%================================================%% 
            %% internal weights updating cycles for critnet   %% 
            %%================================================%% 
                    
            cyc = 0; 
            ecrit = alpha*J-(Jprev-reinf); 
		    Ec = 0.5 * ecrit^2; 
            while (Ec>Tc & cyc<=Ncrit), 
                    gradEcJ=alpha*ecrit; 
    				%----for the first layer(input to hidden layer)----------- 
                    gradqwc1 = [inputs'; newAction]; 
                    for i=1:N_Hidden, 
                        gradJp = wc2(i); 
                        gradpq = 0.5*(1-p(i)^2); 
                        wc1(:,i) = wc1(:,i) - lc*gradEcJ*gradJp*gradpq*gradqwc1; 
                    end        
                    %----for the second layer(hidden layer to output)----------- 
                    gradJwc2=p'; 
                    wc2 = wc2- lc*gradEcJ*gradJwc2; 
                    %----compute new  J---- 
                    inp=[inputs newAction]; 
                    qc=inp*wc1; 
                    p = (1 - exp(-qc))./(1 + exp(-qc)); 
                    J=p*wc2; 
 
                cyc = cyc +1; 
                ecrit = alpha*J-(Jprev-reinf); 
                Ec = 0.5 * ecrit^2; 
            end                                 % end of "while (Ec>0.05 & cyc<=Ncrit)" 
             
            %normalization weights for critical network 
                if (max(max(abs(wc1)))>1.5) 
                    wc1=wc1/max(max(abs(wc1))); 
                end 
                if max(max(abs(wc2)))>1.5 
                    wc2=wc2/max(max(abs(wc2))); 
                end 
                      
            %%=============================================%% 
            %% internal weights updating cycles for actnet %% 
            %%=============================================%%                 
            cyc = 0;             
            eact = J - Uc; 
            Ea = 0.5*eact^2; 
            while (Ea>Ta & cyc<=Nact), 
                graduv = 0.5*(1-newAction^2);              
                gradEaJ = eact; 
                gradJu = 0; 
                   for i=1:N_Hidden, 
                       gradJu = gradJu + wc2(i)*0.5*(1-p(i)^2)*wc1(WC_Inputs,i); 
                   end 
                    %----for the first layer(input to hidden layer)----------- 
                    for (i=1:N_Hidden), 
                        gradvg = wa2(i); 
                        gradgh = 0.5*(1-g(i)^2); 
                        gradhwa1 = inputs'; 
                        wa1(:,i)=wa1(:,i)-la*gradEaJ*gradJu*graduv*gradvg*gradgh*gradhwa1; 
                    end 
                    %----for the second layer(hidden layer to output)----------- 
                    gradvwa2 = g'; 
                    wa2=wa2-la*gradEaJ*gradJu*graduv*gradvwa2; 
                    %----compute new J and newAction-------  
                    ha = inputs*wa1; 
                    g = (1 - exp(-ha))./(1 + exp(-ha)); 
                    va = g*wa2; 
                    newAction = (1 - exp(-va))./(1 + exp(-va)); 
                     
                    inp=[inputs newAction]; 
                    qc=inp*wc1; 
                    p = (1 - exp(-qc))./(1 + exp(-qc)); 
                    J=p*wc2; 
                         
                cyc = cyc+1; 
                eact = J - Uc; 
                Ea = 0.5*eact^2;   
            end                       %end of "while (Ea>Ta & cyc<=Nact)" 
            
            if ~failure 
                Jprev=J; 
            else 
                break;                %another trial 即跳出“while(lfts<Tit),” 
            end 
            lfts=lfts+1; 
        end                           %end of "while(lfts<Tit)" 结束内部循环 
           
        msgstr1=['Trial # ' int2str(trial) ' has  ' int2str(lfts) ' time steps.']; 
   	    msgstr21=['Trial # ' int2str(trial) ' has successfully balanced for at least ']; 
   	    msgstr22=[msgstr21 int2str(lfts) ' time steps ']; 
A_027