1
+ #include " stdafx.h"
2
+ #include " ObjectModel.h"
3
+ #include " NetworkCase.h"
4
+ #include " Exceptions.h"
5
+ #include " Utility.h"
6
+ #include " DCFlowSolver.h"
7
+ #include < algorithm>
8
+ #include < Eigen/Sparse>
9
+ #include < Eigen/SparseLU>
10
+
11
+ using namespace std ;
12
+ using namespace PowerSolutions ::ObjectModel;
13
+ using namespace Eigen ;
14
+
15
+ namespace PowerSolutions
16
+ {
17
+ namespace PowerFlow
18
+ {
19
+ DCFlowSolver::~DCFlowSolver ()
20
+ {
21
+ }
22
+
23
+ shared_ptr<Solution> DCFlowSolver::Solve (NetworkCase& network)
24
+ {
25
+ auto pn = network.ToPrimitive (PrimitiveNetworkOptions::NodeReorder | PrimitiveNetworkOptions::IgnoreShuntAdmittance);
26
+ return this ->Solve (*pn);
27
+ }
28
+
29
+ shared_ptr<Solution> DCFlowSolver::Solve (PrimitiveNetwork& network)
30
+ {
31
+ auto s = make_shared<Solution>();
32
+ s->Status (SolutionStatus::Success);
33
+ s->IterationCount (0 );
34
+ s->MaxDeviation (0 );
35
+ s->NodeCount (network.Nodes ().size ());
36
+ s->PQNodeCount (network.PQNodes ().size ());
37
+ s->PVNodeCount (network.PVNodes ().size ());
38
+ s->SlackNode (network.SlackNode ()->Bus ());
39
+ // 如果在结果生成完毕之前,直接使用 return s;
40
+ // 就会返回一个表示求解失败的结果。
41
+ s->Status (SolutionStatus::IterationFailed);
42
+ // 计算直流潮流
43
+ // P = Bθ
44
+ auto & YMatrix = network.Admittance ;
45
+ // TODO: 考虑此处为何取负。
46
+ auto BMatrix = -YMatrix.imag ();
47
+ // 生成注入功率向量
48
+ auto Nodes = network.Nodes ().size ();
49
+ assert (Nodes > 1 );
50
+ auto IndependentNodes = Nodes - 1 ;
51
+ SparseLU<SparseMatrix<double >> solver;
52
+ // 求解直流潮流。
53
+ VectorXd PVector;
54
+ PVector.resize (IndependentNodes);
55
+ for (auto & node : network.Nodes ())
56
+ {
57
+ if (node->Type () != NodeType::SlackNode) PVector[node->Index ()] = node->ActivePowerInjection ();
58
+ }
59
+ solver.compute (BMatrix.topLeftCorner (IndependentNodes, IndependentNodes));
60
+ if (solver.info () != Eigen::Success) return s;
61
+ // tcout << _T(" B矩阵秩:") << Solver.rank() << endl;
62
+ // 除平衡节点外其它节点的 theta 。
63
+ VectorXd theta = solver.solve (PVector);
64
+ auto ps = make_shared<PrimitiveSolution>(network);
65
+ // cout << ThetaVector << endl;
66
+ unordered_map<Component*, ComponentFlowSolution> componentPowerInjections{};
67
+ // 初始化节点功率
68
+ for (auto & ns : ps->NodeStatus ())
69
+ {
70
+ // 直流潮流,电压幅值为1。
71
+ if (ns.Type () == NodeType::SlackNode)
72
+ ns.SetVoltage (1 , 0 );
73
+ else
74
+ ns.SetVoltage (1 , theta[ns.Index ()]);
75
+ }
76
+ // 计算并写入支路功率。
77
+ for (auto & obj : network.SourceNetwork ()->Objects ())
78
+ {
79
+ auto c = dynamic_cast <Component*>(obj);
80
+ if (c != nullptr )
81
+ {
82
+ // 计算元件的潮流。
83
+ // 注意,此处计算元件潮流时,还是按照交流情况进行计算的。只不过仅保留计算结果的实数部分。
84
+ // 可能会造成结果偏差。
85
+ auto cflow = c->EvalComponentFlow (*ps);
86
+ auto injections = cflow.PowerInjections ();
87
+ for (size_t i = 0 ; i < cflow.PowerInjections ().size (); i++)
88
+ {
89
+ complexd inj (cflow.PowerInjections (i).real (), 0 );
90
+ // 修正潮流结果。
91
+ cflow.PowerInjections (i, inj);
92
+ // 记录潮流结果。
93
+ // 并累加节点注入功率。
94
+ ps->NodeStatus (network.Nodes (c->Buses (i)).Index ()).AddPowerInjections (inj.real (), inj.imag ());
95
+ }
96
+ cflow.PowerShunt (0 );
97
+ // 注意,此处不能直接 AddComponentFlow
98
+ // 因为 AddComponentFlow 会尝试修改 Solution.NodeFlow 的内容。
99
+ // 而此时 NodeFlow 为空。
100
+ componentPowerInjections.emplace (c, move (cflow));
101
+ }
102
+ }
103
+ // 写入节点功率。
104
+ for (auto & node : network.Nodes ()) s->AddNodeFlow (node->Bus (), ps->NodeStatus (node->Index ()));
105
+ // 写入元件注入功率。
106
+ for (auto & p : componentPowerInjections) s->AddComponentFlow (p.first , move (p.second ));
107
+ s->Status (SolutionStatus::Success);
108
+ return s;
109
+ }
110
+ }
111
+ }
0 commit comments