Optimization: Attempting to declare variables in matrix format and exporting in excel format

Question

I have price, cost, demand in different sheets of an excel workbook. Each of these have two dimensions such as Product and Location. Total number of columns(timeperiods) is 6. I am looking to declare my output variables in matrix format so that I am able to export them in excel format. My current variable output is as:

plant__(('Product_1',_'Location_1'),_'Prd_1') 197.0
plant__(('Product_1',_'Location_1'),_'Prd_2') 240.0
plant__(('Product_1',_'Location_1'),_'Prd_3') 167.0
plant__(('Product_2',_'Location_1'),_'Prd_1') 260.0
plant__(('Product_2',_'Location_1'),_'Prd_2') 114.0
plant__(('Product_2',_'Location_1'),_'Prd_3') 253.0

My code is as follows:

from pulp import *
import pandas as pd
price=pd.read_excel("Path\\Python.xlsx", sheet_name='Planned Price', index_col=(0,1))
consen_dem=pd.read_excel("Path\\Python.xlsx", sheet_name='Consensus', index_col=(0,1))
cost=pd.read_excel("Path\\Python.xlsx", sheet_name='Cost', index_col=(0,1))
prod=['Product 1','Product 2','Product 3','Product 4']
time=['Prd 1','Prd 2','Prd 3']
days=pd.read_excel("Path\\Python.xlsx", sheet_name='Days')
product=pd.read_excel("Path\\Python.xlsx", sheet_name='Products', index_col=(0,1))
location=pd.read_excel("Path\\Python.xlsx", sheet_name='Locations')
k = product.shape
lociod=['Location 1','Location 2']    
da = [('Product 1','Location 1'),('Product 2', 'Location 1')]
const=input("Do you want to run constraints or not? (Y/N)")           
#model initialization
m=LpProblem("Model", LpMaximize)
#adding variable               
y = LpVariable.dicts("plant_", [(d,j) for d in da for j in time], lowBound=0, upBound=None, cat='continuous')
#defining objective
m += (lpSum([(price.loc[d,j]-cost.loc[d,j]) * y[(d,j)]   for d in da for j in time])- lpSum([(consen_dem.loc[d,j]-y[(d,j)])*5 for d in da for j in time])) 
#constraints
for d in da:
    for j in time:
        if heur=='Y':
            m+=lpSum(y[(d,j)])>=0
        else:
            m+=lpSum(y[(d,j)])<=200
for d in da:
    for j in time:
        m+=lpSum(y[(d,j)]) <= consen_dem.loc[d][j] 
#execute
m.solve()
#output
dict_prod={}
print('\n' + "Status: {}".format(LpStatus[m.status]))
for i in m.variables():
    dict_prod[i]=i.varValue
        print(i,i.varValue)
print(m.variables())
print("Gross Profit = {:,} ($/Month)".format(int(value(m.objective))))