2-马尔可夫决策过程.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true,
    "ExecuteTime": {
     "end_time": "2023-06-07T10:41:29.775580800Z",
     "start_time": "2023-06-07T10:41:29.541868500Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "outputs": [],
   "source": [
    "np.random.seed(0)\n",
    "# 定义状态转移概率矩阵P\n",
    "P = np.array([\n",
    "    [0.9, 0.1, 0.0, 0.0, 0.0, 0.0],\n",
    "    [0.5, 0.0, 0.5, 0.0, 0.0, 0.0],\n",
    "    [0.0, 0.0, 0.0, 0.6, 0.0, 0.4],\n",
    "    [0.0, 0.0, 0.0, 0.0, 0.3, 0.7],\n",
    "    [0.0, 0.2, 0.3, 0.5, 0.0, 0.0],\n",
    "    [0.0, 0.0, 0.0, 0.0, 0.0, 1.0],\n",
    "])\n",
    "\n",
    "rewards = [-1, -2, -2, 10, 1, 0] # 定义奖励函数R(s)\n",
    "gamma = 0.5 # 定义折扣因子"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-06-07T10:41:33.967055100Z",
     "start_time": "2023-06-07T10:41:33.651957900Z"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "outputs": [],
   "source": [
    "# 给定一条序列，计算从某个索引（起始状态）开始到序列最后（终止状态）得到的回报\n",
    "def compute_return(start_index, chain, gamma):\n",
    "    G = 0\n",
    "    # 用reverse实现 γ^(i-start_index) 是因为贝尔曼方程用于计算状态值或动作值函数的更新，其中的累积奖励（G）的计算需要考虑未来的奖励。\n",
    "    # 反向迭代的方式保证了在计算当前状态或动作的累积奖励时，已经考虑了所有未来奖励的折扣，可以保持一致性和准确性。\n",
    "    for i in reversed(range(start_index, len(chain))):\n",
    "        G = gamma * G + rewards[chain[i] - 1]\n",
    "    return G"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-06-07T10:41:35.233824100Z",
     "start_time": "2023-06-07T10:41:34.989013400Z"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "根据本序列计算得到回报为：-2.5\n"
     ]
    }
   ],
   "source": [
    "# 给一个状态序列，s1-s2-s3-s6\n",
    "chain = [1, 2, 3, 6]\n",
    "start_index = 0\n",
    "G = compute_return(start_index, chain, gamma)\n",
    "print(\"根据本序列计算得到回报为：%s\" % G)"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-06-07T10:41:36.545890800Z",
     "start_time": "2023-06-07T10:41:36.339924700Z"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "outputs": [],
   "source": [
    "# 计算状态价值函数 V = (E - γP)^(-1)*R\n",
    "def compute_value(P, R, gamma, states_num):\n",
    "    R = np.array(R).reshape((-1, 1)) # 展成一列\n",
    "    I = np.eye(states_num) # 生成一个对角矩阵\n",
    "    C = I - np.multiply(gamma, P)\n",
    "    C_1 = np.linalg.inv(C) # 矩阵求逆 用inv和pinv结果不同\n",
    "    return np.dot(C_1, R)"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-06-07T10:41:38.249199400Z",
     "start_time": "2023-06-07T10:41:38.031507400Z"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MRP中每个状态价值为\n",
      " [[-2.01950168]\n",
      " [-2.21451846]\n",
      " [ 1.16142785]\n",
      " [10.53809283]\n",
      " [ 3.58728554]\n",
      " [ 0.        ]]\n"
     ]
    }
   ],
   "source": [
    "V = compute_value(P, rewards, gamma, P.shape[0])\n",
    "print(\"MRP中每个状态价值为\\n\", V)"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-06-07T10:41:39.551349900Z",
     "start_time": "2023-06-07T10:41:39.304964Z"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "outputs": [],
   "source": [
    "# MDP的计算 转化成MRP 用MRP的方法\n",
    "S = [\"s1\", \"s2\", \"s3\", \"s4\", \"s5\"]  # 状态集合\n",
    "A = [\"保持s1\", \"前往s1\", \"前往s2\", \"前往s3\", \"前往s4\", \"前往s5\", \"概率前往\"]  # 动作集合\n",
    "# 状态转移函数p(s'|s,a) 前者是策略决定 后者是环境决定\n",
    "P = {\n",
    "    \"s1-保持s1-s1\": 1.0,\n",
    "    \"s1-前往s2-s2\": 1.0,\n",
    "    \"s2-前往s1-s1\": 1.0,\n",
    "    \"s2-前往s3-s3\": 1.0,\n",
    "    \"s3-前往s4-s4\": 1.0,\n",
    "    \"s3-前往s5-s5\": 1.0,\n",
    "    \"s4-前往s5-s5\": 1.0,\n",
    "    \"s4-概率前往-s2\": 0.2,\n",
    "    \"s4-概率前往-s3\": 0.4,\n",
    "    \"s4-概率前往-s4\": 0.4,\n",
    "}\n",
    "# 奖励函数 这里定义成r(s,a) 这意味着无论从状态s转移到哪个状态s' 奖励的大小都是相同的\n",
    "R = {\n",
    "    \"s1-保持s1\": -1,\n",
    "    \"s1-前往s2\": 0,\n",
    "    \"s2-前往s1\": -1,\n",
    "    \"s2-前往s3\": -2,\n",
    "    \"s3-前往s4\": -2,\n",
    "    \"s3-前往s5\": 0,\n",
    "    \"s4-前往s5\": 10,\n",
    "    \"s4-概率前往\": 1,\n",
    "}\n",
    "gamma = 0.5  # 折扣因子\n",
    "MDP = (S, A, P, R, gamma)\n",
    "\n",
    "# 策略1,随机策略\n",
    "Pi_1 = {\n",
    "    \"s1-保持s1\": 0.5,\n",
    "    \"s1-前往s2\": 0.5,\n",
    "    \"s2-前往s1\": 0.5,\n",
    "    \"s2-前往s3\": 0.5,\n",
    "    \"s3-前往s4\": 0.5,\n",
    "    \"s3-前往s5\": 0.5,\n",
    "    \"s4-前往s5\": 0.5,\n",
    "    \"s4-概率前往\": 0.5,\n",
    "}\n",
    "# 策略2\n",
    "Pi_2 = {\n",
    "    \"s1-保持s1\": 0.6,\n",
    "    \"s1-前往s2\": 0.4,\n",
    "    \"s2-前往s1\": 0.3,\n",
    "    \"s2-前往s3\": 0.7,\n",
    "    \"s3-前往s4\": 0.5,\n",
    "    \"s3-前往s5\": 0.5,\n",
    "    \"s4-前往s5\": 0.1,\n",
    "    \"s4-概率前往\": 0.9,\n",
    "}\n",
    "\n",
    "\n",
    "# 把输入的两个字符串通过\"-\"连接 便于使用上述定义的P、R变量\n",
    "def join(str1, str2):\n",
    "    return str1 + '-' + str2"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-06-07T10:41:41.592164600Z",
     "start_time": "2023-06-07T10:41:41.342768900Z"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MRP中每个状态价值为\n",
      " [[-1.22555411]\n",
      " [-1.67666232]\n",
      " [ 0.51890482]\n",
      " [ 6.0756193 ]\n",
      " [ 0.        ]]\n"
     ]
    }
   ],
   "source": [
    "# 本例中给定的是R(s',s,a) 故由Pi_1和P的MDP转化得的MRP 用MRP的方式进行计算\n",
    "MRP = [\n",
    "    [0.5, 0.5, 0.0, 0.0, 0.0],\n",
    "    [0.5, 0.0, 0.5, 0.0, 0.0],\n",
    "    [0.0, 0.0, 0.0, 0.5, 0.5],\n",
    "    [0.0, 0.1, 0.2, 0.2, 0.5],\n",
    "    [0.0, 0.0, 0.0, 0.0, 1.0],\n",
    "]\n",
    "MRP = np.array(MRP)\n",
    "# 由R和Pi_1得到的奖励矩阵\n",
    "R = [-0.5, -1.5, -1.0, 5.5, 0]\n",
    "\n",
    "V = compute_value(MRP, R, gamma, len(S))\n",
    "print(\"MRP中每个状态价值为\\n\", V)"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-06-07T10:41:43.243799800Z",
     "start_time": "2023-06-07T10:41:42.988051Z"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "outputs": [],
   "source": [
    "def sample(MDP, Pi, timestep_max, number):\n",
    "    \"\"\"蒙特卡洛法采样序列\n",
    "        MDP 马尔可夫参数\n",
    "        Pi 策略\n",
    "        timestep_max 限制最长时间步\n",
    "        number 总共采样序列数\n",
    "    \"\"\"\n",
    "    S, A, P, R, gamma = MDP\n",
    "    episodes = []\n",
    "    # 采样number次\n",
    "    for _ in range(number):\n",
    "        episode = []\n",
    "        timestep = 0\n",
    "        # 随机初始化一个状态s0 这里不让它一开始就是s5\n",
    "        s = S[np.random.randint(0, 4)]\n",
    "        # 如果采样数达到最大或到达最终状态s5 停止\n",
    "        while timestep < timestep_max and s != \"s5\":\n",
    "            timestep += 1\n",
    "            rand, temp = np.random.rand(), 0\n",
    "            # 由策略&当前状态 选择当前动作\n",
    "            for a_opt in A:\n",
    "                temp += Pi.get(join(s,a_opt), 0)\n",
    "                if rand < temp:\n",
    "                    a = a_opt\n",
    "                    r = R.get(join(s, a_opt), 0)\n",
    "                    break\n",
    "            # 由状态转移矩阵&当前状态&当前动作 选择下一个状态\n",
    "            rand, temp = np.random.rand(), 0\n",
    "            for s_opt in S:\n",
    "                temp += P.get(join(join(s,a), s_opt), 0)\n",
    "                if rand < temp:\n",
    "                    s_next = s_opt\n",
    "                    break\n",
    "            # 放到池子里\n",
    "            episode.append((s, a, r, s_next))\n",
    "            s = s_next\n",
    "        episodes.append(episode)\n",
    "    return episodes"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-06-07T10:41:45.322619500Z",
     "start_time": "2023-06-07T10:41:45.099746300Z"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "第一条序列：\n",
      " [('s1', '前往s2', 0, 's2'), ('s2', '前往s3', -2, 's3'), ('s3', '前往s5', 0, 's5')]\n",
      "第二条序列：\n",
      " [('s4', '概率前往', 1, 's4'), ('s4', '前往s5', 10, 's5')]\n",
      "第三条序列：\n",
      " [('s4', '前往s5', 10, 's5')]\n",
      "第四条序列：\n",
      " [('s2', '前往s1', -1, 's1'), ('s1', '保持s1', -1, 's1'), ('s1', '前往s2', 0, 's2'), ('s2', '前往s3', -2, 's3'), ('s3', '前往s4', -2, 's4'), ('s4', '前往s5', 10, 's5')]\n",
      "第五条序列：\n",
      " [('s2', '前往s3', -2, 's3'), ('s3', '前往s4', -2, 's4'), ('s4', '前往s5', 10, 's5')]\n"
     ]
    }
   ],
   "source": [
    "# 采样5次 最大20步 示例\n",
    "timestep_max, number = 20, 5\n",
    "episodes = sample(MDP, Pi_1, timestep_max, number)\n",
    "print(\"第一条序列：\\n\", episodes[0])\n",
    "print(\"第二条序列：\\n\", episodes[1])\n",
    "print(\"第三条序列：\\n\", episodes[2])\n",
    "print(\"第四条序列：\\n\", episodes[3])\n",
    "print(\"第五条序列：\\n\", episodes[4])"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-06-07T10:41:47.547657300Z",
     "start_time": "2023-06-07T10:41:47.281309100Z"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "outputs": [],
   "source": [
    "def MC(episodes, V, N, gamma):\n",
    "    \"\"\"\n",
    "        蒙特卡洛法，每当s_i出现，计算对应的G_t，累加后取均值\n",
    "    \"\"\"\n",
    "    for episode in episodes:\n",
    "        G = 0\n",
    "        # 从后往前计算\n",
    "        for i in range(len(episode) - 1, -1, -1):\n",
    "            (s, a, r, s_next) = episode[i]\n",
    "            G = r + gamma * G\n",
    "            N[s] += 1\n",
    "            # 和对比老虎机一样 常用的累加计算方式\n",
    "            V[s] += (G - V[s]) / N[s]"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-06-07T10:42:36.832553100Z",
     "start_time": "2023-06-07T10:42:36.601181Z"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "使用蒙特卡洛方法计算MDP的状态价值为\n",
      " {'s1': -1.213994389377003, 's2': -1.6519564617977478, 's3': 0.5647939497667427, 's4': 6.193547684776152, 's5': 0}\n"
     ]
    }
   ],
   "source": [
    "# 采样1000次 每次20步 可以自行修改\n",
    "episodes = sample(MDP, Pi_1, 20, 1000)\n",
    "gamma = 0.5\n",
    "V = {\"s1\": 0, \"s2\": 0, \"s3\": 0, \"s4\": 0, \"s5\": 0}\n",
    "N = {\"s1\": 0, \"s2\": 0, \"s3\": 0, \"s4\": 0, \"s5\": 0}\n",
    "MC(episodes, V, N, gamma)\n",
    "print(\"使用蒙特卡洛方法计算MDP的状态价值为\\n\", V)"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-06-07T10:42:38.288306100Z",
     "start_time": "2023-06-07T10:42:37.914613600Z"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "import copy\n",
    "\n",
    "\n",
    "class CliffWalkingEnv:\n",
    "    \"\"\" 悬崖漫步环境\"\"\"\n",
    "    def __init__(self, ncol=12, nrow=4):\n",
    "        self.ncol = ncol  # 定义网格世界的列\n",
    "        self.nrow = nrow  # 定义网格世界的行\n",
    "        # 转移矩阵P[state][action] = [(p, next_state, reward, done)]包含下一个状态和奖励\n",
    "        self.P = self.createP()\n",
    "\n",
    "    def createP(self):\n",
    "        # 初始化\n",
    "        P = [[[] for j in range(4)] for i in range(self.nrow * self.ncol)]\n",
    "        # 4种动作 change[0]:上 change[1]:下 change[2]:左 change[3]:右 坐标系原点(0,0)\n",
    "        # 定义在左上角\n",
    "        change = [[0, -1], [0, 1], [-1, 0], [1, 0]]\n",
    "        for i in range(self.nrow):\n",
    "            for j in range(self.ncol):\n",
    "                for a in range(4):\n",
    "                    # 位置在悬崖或者目标状态因为无法继续交互 任何动作奖励都为0\n",
    "                    if i == self.nrow - 1 and j > 0:\n",
    "                        P[i * self.ncol + j][a] = [(1, i * self.ncol + j, 0,\n",
    "                                                    True)]\n",
    "                        continue\n",
    "                    # 其他位置\n",
    "                    next_x = min(self.ncol - 1, max(0, j + change[a][0]))\n",
    "                    next_y = min(self.nrow - 1, max(0, i + change[a][1]))\n",
    "                    next_state = next_y * self.ncol + next_x\n",
    "                    reward = -1\n",
    "                    done = False\n",
    "                    # 下一个位置在悬崖或者终点\n",
    "                    if next_y == self.nrow - 1 and next_x > 0:\n",
    "                        done = True\n",
    "                        if next_x != self.ncol - 1:  # 下一个位置在悬崖\n",
    "                            reward = -100\n",
    "                    P[i * self.ncol + j][a] = [(1, next_state, reward, done)]\n",
    "        return P"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "class PolicyIteration:\n",
    "    \"\"\" 策略迭代算法 \"\"\"\n",
    "    def __init__(self, env, theta, gamma):\n",
    "        self.env = env\n",
    "        self.v = [0] * self.env.ncol * self.env.nrow  # 初始化价值为0\n",
    "        self.pi = [[0.25, 0.25, 0.25, 0.25]\n",
    "                   for i in range(self.env.ncol * self.env.nrow)]  # 初始化为均匀随机策略\n",
    "        self.theta = theta  # 策略评估收敛阈值\n",
    "        self.gamma = gamma  # 折扣因子\n",
    "\n",
    "    def policy_evaluation(self):  # 策略评估\n",
    "        cnt = 1  # 计数器\n",
    "        while 1:\n",
    "            max_diff = 0\n",
    "            new_v = [0] * self.env.ncol * self.env.nrow\n",
    "            for s in range(self.env.ncol * self.env.nrow):\n",
    "                qsa_list = []  # 开始计算状态s下的所有Q(s,a)价值\n",
    "                for a in range(4):\n",
    "                    qsa = 0\n",
    "                    for res in self.env.P[s][a]:\n",
    "                        p, next_state, r, done = res\n",
    "                        qsa += p * (r + self.gamma * self.v[next_state] * (1 - done))\n",
    "                    qsa_list.append(self.pi[s][a] * qsa)\n",
    "                new_v[s] = sum(qsa_list)  # 状态价值函数和动作价值函数之间的关系\n",
    "                max_diff = max(max_diff, abs(new_v[s] - self.v[s]))\n",
    "            self.v = new_v\n",
    "            if max_diff < self.theta: break  # 满足收敛条件 退出评估迭代\n",
    "            cnt += 1\n",
    "        print(\"策略评估进行%d轮后完成\" % cnt)\n",
    "\n",
    "    def policy_improvement(self):  # 策略提升\n",
    "        for s in range(self.env.nrow * self.env.ncol):\n",
    "            qsa_list = []\n",
    "            for a in range(4):\n",
    "                qsa = 0\n",
    "                for res in self.env.P[s][a]:\n",
    "                    p, next_state, r, done = res\n",
    "                    qsa += p * (r + self.gamma * self.v[next_state] * (1 - done))\n",
    "                qsa_list.append(qsa)\n",
    "            maxq = max(qsa_list)\n",
    "            cntq = qsa_list.count(maxq)  # 计算有几个动作得到了最大的Q值\n",
    "            # 让这些动作均分概率\n",
    "            self.pi[s] = [1 / cntq if q == maxq else 0 for q in qsa_list]\n",
    "        print(\"策略提升完成\")\n",
    "        return self.pi\n",
    "\n",
    "    def policy_iteration(self):  # 策略迭代\n",
    "        while 1:\n",
    "            self.policy_evaluation()\n",
    "            old_pi = copy.deepcopy(self.pi)  # 将列表进行深拷贝,方便接下来进行比较\n",
    "            new_pi = self.policy_improvement()\n",
    "            if old_pi == new_pi: break"
   ],
   "metadata": {
    "collapsed": false
   }
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}