-
Notifications
You must be signed in to change notification settings - Fork 3
/
sens_chart.py
128 lines (111 loc) · 3.88 KB
/
sens_chart.py
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
" create bar chart "
import sys
import matplotlib.pyplot as plt
import numpy as np
from numpy import random
from gpkit.repr_conventions import unitstr
from solar import Mission, Aircraft
#pylint: disable=invalid-name, anomalous-backslash-in-string
def get_highestsens(model, res, varnames=None, N=10):
" plot bar chart of sensitivities "
pss = []
ngs = []
sens = {}
if varnames:
for vname in varnames:
sen = res["sensitivities"]["constants"][vname]
if hasattr(sen, "__len__"):
print vname
print sen
val = max(np.abs(sen.values()))
vk = [svk for svk in sen if abs(sen[svk]) == val][0]
sen = sum(sen.values())
else:
vk = model[vname].key
sens[vk] = sen
else:
for s in res["sensitivities"]["constants"]:
val = res["sensitivities"]["constants"][s]
if hasattr(val, "shape"):
if len(val.shape) == 1:
val = sum(val)
if len(val.shape) == 2:
val = sum(np.hstack(val))
sens[model[s].key] = val
labels = []
i = 0
sorted_sens = dict_sort(sens)
for s in sorted_sens:
if i > N:
break
i += 1
vk = s[0]
val = sum(np.hstack([res(vk)]))
if "units" in vk.descr:
uts = unitstr(vk.descr["units"])
else:
uts = ""
lbl = vk.descr["label"]
labels.append(lbl + "$ =%.2f$ %s" % (val, uts.replace("*", "")))
if s[1] > 0:
pss.append(s[1])
ngs.append(0)
else:
ngs.append(abs(s[1]))
pss.append(0)
ind = np.arange(0.5, i + 0.5, 1)
sensdict = {"positives": pss, "negatives": ngs, "indicies": ind,
"labels": labels}
return sensdict
def dict_sort(vdict):
" sort variable sensitivity dict"
slist = [(0, 0.0)]
for v in vdict:
for i, sv in enumerate(slist):
if abs(vdict[v]) > abs(sv[1]):
slist.insert(i, (v, vdict[v]))
break
del slist[-1]
return slist
def plot_chart(sensdict):
"plot sensitivities on bar chart"
fig, ax = plt.subplots()
ax.bar(sensdict["indicies"], sensdict["positives"], 0.5, color="#4D606E")
ax.bar(sensdict["indicies"], -np.array(sensdict["negatives"]), 0.5,
color="#3FBAC2")
ax.set_xlim([0.0, sensdict["indicies"][-1]+1])
ax.set_xticks(sensdict["indicies"])
ax.set_xticklabels(sensdict["labels"], rotation=-45, ha="left")
# ax.legend(["Positive", "Negative"])
ax.set_ylabel("sensitivities")
ax.grid()
return fig, ax
def test():
" test for integrated testing "
v = Aircraft(Npod=0, sp=True)
model = Mission(v, latitude=[20])
model.cost = model[model.aircraft.Wtotal]
result = model.localsolve("mosek")
_ = get_highestsens(model, result)
vn = {model.aircraft.Wpay: "$W_{\\mathrm{pay}}$",
model.aircraft.battery.etacharge: "$\\eta_{\\mathrm{charge}}$"}
_ = get_highestsens(model, result, vn)
if __name__ == "__main__":
if len(sys.argv) > 1:
path = sys.argv[1]
else:
path = ""
V = Aircraft(Npod=0, sp=False)
M = Mission(V, latitude=[20])
M.cost = M[M.aircraft.Wtotal]
sol = M.solve("mosek")
vns = {M.aircraft.Wpay: "$W_{\\mathrm{pay}}$",
M.aircraft.battery.etacharge: "$\\eta_{\\mathrm{charge}}$",
M.aircraft.battery.etadischarge: "$\\eta_{\\mathrm{discharge}}$",
M.aircraft.battery.hbatt: "$h_{\\mathrm{batt}}$",
M.aircraft.solarcells.etasolar: "$\\eta_{\\mathrm{solar}}$",
M.mission[1].winggust.Nmax: "$N_{\\mathrm{load}}$",
M.mission[1].aircraftPerf.drag.wing.e: "$e$"}
sd = get_highestsens(M, sol, vns)
f, a = plot_chart(sd)
f.savefig(path + "sensbar.pdf", bbox_inches="tight")