StoDyM_MOEA.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import time\n",
    "from pandas import read_excel, DataFrame, ExcelWriter\n",
    "from collections import OrderedDict\n",
    "from datetime import timedelta\n",
    "from stodymmoea import *\n",
    "from platypus import *\n",
    "\n",
    "\n",
    "def pareto_front_dataframe(pf_l):\n",
    "    \n",
    "    # empty Pareto-front dataframe\n",
    "    pf_df = DataFrame(columns=pfdfcn_l)\n",
    "    \n",
    "    # Pareto-front's dam targeted location list\n",
    "    pfdtl_l = [[dv_l[0].decode(v1.variables[v2]) for v2 in range(d_no)] for v1 in pf_l]\n",
    "    \n",
    "    # Pareto-front's dam storage capacity list\n",
    "    if dscdv_i == True:  # if Integer decision variable\n",
    "        pfdsc_l = [[dsc_mult*dv_l[1].decode(v1.variables[v2]) for v2 in range(d_no, dv_no)] for v1 in pf_l]\n",
    "    else:  # if Subset decision variable\n",
    "        pfdsc_l = [v1.variables[-1] for v1 in pf_l]\n",
    "    \n",
    "    # Pareto-front's objective value list\n",
    "    pov_l = np.array([[v1.objectives[v2] for v2 in range(o_no)] for v1 in pf_l])\n",
    "    \n",
    "    # Store Pareto-front in dataframe\n",
    "    for i_v1, v1 in enumerate(pfdtl_l): \n",
    "        dtl_a, dtlri_a = np.unique(pfdtl_l[i_v1], return_index=True)\n",
    "        disc_a = np.array(pfdsc_l[i_v1])[dtlri_a]\n",
    "        dda_a = dam_drainage_area(fa_a, dtl_a, sluc_a, slco_a, pa_v)\n",
    "        dda_tv = 100 * (np.sum(dda_a)/ca_t)\n",
    "        dnsv_a, dle_a, ssdr_a = dam_deployment_output(sluc_a, slco_a, fa_a, ssi_a, \n",
    "                                                      dtl_a, dda_a, disc_a, pa_v, ste_tv, dle_max, csw_kg, s_sg)\n",
    "        pddsc_l = [*dtl_a, *disc_a, np.sum(disc_a), dda_tv, *pov_l[i_v1], \n",
    "                   *dle_a, np.mean(dle_a), \n",
    "                   (np.sum(np.power(pov_l[i_v1]-u_p, 2)))**0.5]\n",
    "        pf_df.loc[i_v1, pfdfcn_l] = pddsc_l\n",
    "    \n",
    "    # sorting Pareto-front by dams' targeted location and storage capacity\n",
    "    pf_df = pf_df.drop_duplicates().sort_values(by=didcn_l+disccn_l).reset_index(drop=True) \n",
    "    \n",
    "    return pf_df\n",
    "\n",
    "\n",
    "def stream_connection_structure():\n",
    "    \n",
    "    \"\"\"\n",
    "    This function calculates: \n",
    "    (a) all the connnected upstream segments for each stream segment,\n",
    "    (b) connected order of stream segments from downstream to upstream\n",
    "    \"\"\"\n",
    "    \n",
    "    # an empty array to store stream link upstream connection\n",
    "    sluc_a = np.empty(len(sid_a), dtype=object)\n",
    "    \n",
    "    # an empty list to store connected order of stream segments from downstream to upstream\n",
    "    slco_l = []  \n",
    "    \n",
    "    # end segment\n",
    "    sl_mds = sid_a[np.where(nsid_a == (len(nsid_a) + 1))[0]]\n",
    "    sl_mdsi = np.asscalar(np.where(sid_a == sl_mds)[0])\n",
    "    sluc_a[sl_mdsi] = sl_mds\n",
    "    slco_l.append(sl_mds)\n",
    "    \n",
    "    while True:  \n",
    "        if len(sl_mds) == 1:  # detemine next downstream segments of end segment and store the values\n",
    "            sl_nuci = np.where(nsid_a == sl_mds)[0]\n",
    "            sl_mds = sid_a[sl_nuci]\n",
    "            sluc_a[sl_mdsi] = np.append(sluc_a[sl_mdsi], sl_mds)\n",
    "            slco_l.append(sl_mds) \n",
    "        else:  # update next donwstream segments and store the values\n",
    "            sl_nuci = np.array([np.where(nsid_a == v1)[0] for v1 in sl_mds])\n",
    "            sl_nuc = [sid_a[v2] if len(v2) != 0 else v2 for v2 in sl_nuci] \n",
    "            for v3, v4 in zip(sl_mds, sl_nuc):\n",
    "                if len(v4) != 0:\n",
    "                    sluc_a[v3-1] = np.append(sid_a[v3-1], v4)\n",
    "                else:\n",
    "                    sluc_a[v3-1] = v4\n",
    "            for v5 in sl_mds:  # iterate updated next downstream segments and store the values\n",
    "                sl_uci = np.array([i_v6 for i_v6, v6 in enumerate(sluc_a) \n",
    "                                   if v6 is not None and v5 in v6 and i_v6 != v5 - 1])\n",
    "                for v7 in sl_uci:\n",
    "                    v7_uuc = np.append(sluc_a[v7], sluc_a[v5-1])\n",
    "                    v7_usp = np.unique(v7_uuc, return_index=True)\n",
    "                    sluc_a[v7] = v7_uuc[np.sort(v7_usp[1])] \n",
    "            sl_nenuc = [sid_a[v8] for v8 in sl_nuci if len(v8) != 0]\n",
    "            if len(sl_nenuc) != 0:  # continue if non-empty next upstream connection is found\n",
    "                sl_mds = np.concatenate(sl_nenuc)\n",
    "                slco_l.append(sl_mds)\n",
    "            else:  # otherwise break the loop \n",
    "                break \n",
    "    \n",
    "    # connected order of stream segments from downstream to upstream\n",
    "    slco_a = np.array(slco_l)[::-1]  \n",
    "    \n",
    "    return sluc_a, slco_a  # return the outputs\n",
    "\n",
    "\n",
    "pa_v = (20**2)/10**6  # pixel area value in km^2\n",
    "ste_tv = 0.1  # Sediment Trapping Efficiency threshold value\n",
    "dle_min, dle_max = (20, 30)  # dam life expectancy minimum and maximum values\n",
    "dsc_mult = 10**3  # dam initial storage capacity multiplier for Integer decision variables\n",
    "csw_kg, s_sg = (907.185, 2650)  # connversion of sediment unit to Kg and sediment specific graity\n",
    "eps_v = 0.001  # epsilon value\n",
    "\n",
    "cpn_d = {\"Majiagou\": 34915, \n",
    "         \"Shejiagou\": 12419}  # catchment pixel number dictionary\n",
    "cdn_d = {\"Majiagou\": 10, \n",
    "         \"Shejiagou\": 5}  # catchment dam number dictionary\n",
    "cdscul_d = {\"Majiagou\": 700000, \n",
    "            \"Shejiagou\": 240000}  # catchment dam storage capacity upper limit dictionary\n",
    "cdsci_d = {\"Majiagou\": Integer(40, 100), \n",
    "           \"Shejiagou\": Integer(20, 60)}  # catchment dam storage capacity Integer dictionary\n",
    "cdscs_d = {\"Majiagou\": list(range(40000, 100000+1, 5000)), \n",
    "           \"Shejiagou\": list(range(20000, 60000+1, 1000))}  # catchment dam storage capacity Subset dictionary\n",
    "\n",
    "\n",
    "cn_s = \"Shejiagou\"  # catchment string name\n",
    "ca_t = cpn_d[cn_s] * pa_v  # catchment area (km^2)\n",
    "d_no = cdn_d[cn_s]  # number of dams\n",
    "dsc_ul = cdscul_d[cn_s]  # dam initial storage capacity upper limit\n",
    "\n",
    "# stream characteristics \n",
    "sc_fn = cn_s + \"_Stream_Characteristics.xlsx\"  # stream characteristics file name\n",
    "sc_df = read_excel(sc_fn, sheet_name=\"SC_Data\")  # stream characteristics dataframe\n",
    "sid_a = np.array(sc_df[\"SID\"], dtype=int)  # stream ID array\n",
    "nsid_a = np.array(sc_df[\"NSID\"], dtype=int)  # next stream ID array\n",
    "fa_a = np.array(sc_df[\"DFA\"], dtype=int)  # flow accumulation array\n",
    "ssi_a = np.array(sc_df[\"SSI\"], dtype=float)  # stream sediment input array\n",
    "\n",
    "# stream link upstream connection and connected order\n",
    "sluc_a, slco_a = stream_connection_structure()  \n",
    "\n",
    "# objective, population size, function evaluations, seeds, computer core\n",
    "o_no, ps_no, fe_no, s_no, cc_no = (4, 50, 1000, 4, 4)  \n",
    "\n",
    "dscdv_i = True  # Integer type of dam storage capacity is True\n",
    "\n",
    "# decision variables and other feature \n",
    "if dscdv_i == True:  # if dams' storage capacity is Integer\n",
    "    dv_no = 2 * d_no  # decision variables number\n",
    "    dv_l = [Integer(1, len(sid_a)), cdsci_d[cn_s]]  #  decision variable list\n",
    "    epsnsga2_d = {\"population_size\": ps_no, \n",
    "                  \"epsilons\": eps_v}  # EpsNSGAII features dictionary\n",
    "else:  # if dams' storage capacity is Subset\n",
    "    dv_no = d_no + 1  # decision variables number\n",
    "    dv_l = [Integer(1, len(sid_a)), Subset(cdscs_d[cn_s], d_no)]  #  decision variable list\n",
    "    epsnsga2_d = {\"population_size\": ps_no, \n",
    "                  \"epsilons\": eps_v, \n",
    "                  \"variator\": CompoundOperator(HUX(), BitFlip(), SSX(), Replace())}  # EpsNSGAII features dictionary\n",
    "    \n",
    "od_l = [Problem.MAXIMIZE, Problem.MINIMIZE, Problem.MINIMIZE, Problem.MINIMIZE]  # objective decision list\n",
    "cs_l = [\"==0\", \"==0\", \"==0\", \"<=0\"]  # constraint sign list\n",
    "u_p = np.array([1, 0, 0, 0])  # Utopian point\n",
    "\n",
    "# optimization required values and array\n",
    "ov_t = (dv_l, od_l, cs_l, d_no, dle_min, dle_max, dsc_ul, dsc_mult, pa_v, ste_tv, csw_kg, s_sg)\n",
    "oa_t = (sluc_a, slco_a, fa_a, ssi_a)\n",
    "\n",
    "# optimization output and time\n",
    "if __name__ == \"__main__\":\n",
    "    ss_t = time.time()  \n",
    "    oa_l = [(EpsNSGAII, epsnsga2_d)]  \n",
    "    op_l = [StoDyM(dv_no, o_no, len(cs_l), ov_t, oa_t)]  \n",
    "    with ProcessPoolEvaluator(cc_no) as mp_c:  \n",
    "        or_d = experiment(oa_l, op_l, nfe=fe_no, seeds=s_no, evaluator=mp_c, display_stats=True)  \n",
    "        se_t = time.time()  \n",
    "        r_t = se_t - ss_t\n",
    "\n",
    "# column name list\n",
    "didcn_l = [\"D\"+str(v1) for v1 in range(d_no)]  # dam stream ID\n",
    "disccn_l = [\"SC\"+str(v1)+\"(m^3)\" for v1 in range(d_no)]  # dam storage capacity\n",
    "ocn_l = [\"O-LE\", \"O-SD\", \"O-SDR,ST\", \"O-SDR,LT\"]  # objective\n",
    "dlecn_l = [\"LE\"+str(v1)+\"(yr)\" for v1 in range(d_no)]  # dam life expectancy\n",
    "pfdfcn_l = [*didcn_l, *disccn_l, \"ITSC(m^3)\", \"ITDA(%)\", *ocn_l,  \n",
    "            *dlecn_l, \"LE_A(yr)\", \"Distance\"]  # Pareto-front dataframe\n",
    "\n",
    "# write the output\n",
    "for v1 in or_d.keys():  # iterate algorithm list\n",
    "    for v2 in or_d[v1].keys():  # iterate problem list\n",
    "        v2_fn = \"_\".join(map(str, [cn_s, v1, ps_no, \n",
    "                                   str(int(fe_no/10**3))+\"k\", \"LS\"]))  # simulation output file name\n",
    "        v2xl_w = ExcelWriter(v2_fn+\".xlsx\", engine=\"xlsxwriter\")  # Excel writer object\n",
    "        v2ms_l = []  # an empty list to merge the output of all seeds\n",
    "        for i_v3, v3 in enumerate(or_d[v1][v2]):  # iterate seed solution \n",
    "            v3fs_l = [v4 for v4 in v3 if v4.feasible]  # feasible solution list\n",
    "            v2ms_l.extend(v3fs_l)  # merge feasible solutions\n",
    "        v2pf_l = nondominated(v2ms_l)  # list of non-dominated solutions\n",
    "        v2pf_df = pareto_front_dataframe(v2pf_l)  # Pareto-front dataframe with detailed charateristics\n",
    "        v2pf_df.to_excel(v2xl_w, sheet_name=\"Pareto_Front\", index=False)  # write Pareto-efficient set \n",
    "        t_df = DataFrame(OrderedDict([(\"Total\", r_t), \n",
    "                                      (\"HMS_Format\", str(timedelta(seconds=round(r_t))))]), \n",
    "                         index=[0]).T  # time dataframe\n",
    "        t_df.to_excel(v2xl_w, sheet_name=\"Time\", index=True, header=False)  # write time\n",
    "        v2xl_w.save()  # save the excel file"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}