fengqi/fq_sim
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

add a new method for NR loop, only diode have finished,

Browse Source
This commit is contained in:
feng-arch 2024-08-22 11:37:28 +08:00
parent 633fd2aa7e
commit 10bca9eafc
4 changed files with 177 additions and 76 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
inc/Simulator.hpp
View File

@ -16,6 +16,8 @@ public:
int m;//MNA矩阵比NA矩阵多的行数
Simulator(Circuit c);
Vector dc_solve();
// 上面这个dc_solve()函数存在问题，虽然可以求出结果，但是不是NR迭代
Vector dc_solve_new();
int time_solve( double dt, double tmax, std::string filename);
void add_lc(double dt);
void load_static_G();
@ -32,4 +34,5 @@ void test_mosfet();
void test_vccs();
void test_l6b();
void test_vcr();
void test_dcsolver_new();
#endif // !

2
main.cpp
View File

@ -6,6 +6,6 @@
// Example usage
int main()
{
test_vcr();
test_dcsolver_new();
return 0;
}

242
src/Simulator.cpp
View File

@ -11,60 +11,44 @@ Vector Simulator::dc_solve()
// 牛顿-拉夫森迭代法求解非线性电路
int iter_nums = 0;
double err = 1;
while (err > 1e-9)
while (err > 1e-8)
{
Matrix G_iter = G; // 复制导纳矩阵
Vector I_iter = I; // 复制电流矢量
// 处理vcr器件
// int vcrI_tempndex = n + circuit.vcvs.size() + circuit.cccs.size() + circuit.ccvs.size() * 2 + circuit.voltageSources.size() + circuit.inductors.size(); // 从电感之后开始
int vcrI_tempndex = n + circuit.vcvs.size() + circuit.cccs.size() + circuit.ccvs.size() * 2 + circuit.voltageSources.size() + circuit.inductors.size(); // 从电感之后开始
for (const auto &vcr : circuit.vcrs)
{
int cp = vcr.controlNodePos - 1;
int cn = vcr.controlNodeNeg - 1;
int op = vcr.outputNodePos - 1;
int on = vcr.outputNodeNeg - 1;
// G=1/kU
double U = (cp >= 0 ? V[cp] : 0) - (cn >= 0 ? V[cn] : 0);
double G = get_vcr_G(vcr, U);
// double out_U = 0;
// out_U = (op >= 0 ? V[op] : 0) - (on >= 0 ? V[on] : 0);
// if(G==0)
// {
// out_U = 0;
// }
// else
// {
// out_U = -V[vcrI_tempndex]/G;
// }
// if (op >= 0 && cp >= 0)
// {
// G_iter[op][cp] += G;
// }
// if (op >= 0 && cn >= 0)
// {
// G_iter[op][cn] -= G;
// }
// if (on >= 0 && cp >= 0)
// {
// G_iter[on][cp] -= G;
// }
// if (on >= 0 && cn >= 0)
// {
// G_iter[on][cn] += G;
// }
if(op>=0)
double control_U = (cp >= 0 ? V[cp] : 0) - (cn >= 0 ? V[cn] : 0);
// double out_U = -(op >= 0 ? V[op] : 0) - (on >= 0 ? V[on] : 0);
if (cp >= 0)
{
G_iter[op][op] += G;
G_iter[cp][vcrI_tempndex] -= vcr.transresistance * V[vcrI_tempndex];
}
if(on>=0)
if (cn >= 0)
{
G_iter[on][on] += G;
}
if(op>0&&on>0)
{
G_iter[op][on] -= G;
G_iter[on][op] -= G;
G_iter[cn][vcrI_tempndex] += vcr.transresistance * V[vcrI_tempndex];
}
G_iter[vcrI_tempndex][vcrI_tempndex] -= vcr.transresistance * control_U;
I_iter[vcrI_tempndex] += V[vcrI_tempndex] * vcr.transresistance * control_U;
// 下面的方式使用电阻的方式stamp，但是效果不是很理想
// if (op >= 0)
// {
// G_iter[op][op] += G;
// }
// if (on >= 0)
// {
// G_iter[on][on] += G;
// }
// if (op > 0 && on > 0)
// {
// G_iter[op][on] -= G;
// G_iter[on][op] -= G;
// }
// I_iter[vcrI_tempndex] = out_U;
// vcrI_tempndex++;
}
@ -177,11 +161,14 @@ Vector Simulator::dc_solve()
V = 0.1 * (V_new - V) + V;
continue;
}
else
{
V = 0.9 * V_new + 0.1 * V;
}
if (iter_nums > 1000)
{
break;
}
V = V_new;
}
Vector result(n + m + 2);
result[0] = iter_nums;
@ -267,8 +254,9 @@ void Simulator::load_static_G()
m += circuit.ccvs.size() * 2 + circuit.cccs.size();
m += circuit.voltageSources.size();
m += circuit.inductors.size();
// m += circuit.vcrs.size();
m += circuit.vcrs.size();
Matrix G_temp(n + m, n + m); // 导纳矩阵（扩大以包含辅助变量)
// std::cout << n << " " << m << std::endl;
// 填充导纳矩阵G_temp和电流矢量I_temp
for (const auto &resistor : circuit.resistors)
{
@ -420,25 +408,25 @@ void Simulator::load_static_G()
dcI_tempndex++;
}
// 添加VCRS的影响
// int vcrI_tempndex = n + circuit.vcvs.size() + circuit.cccs.size() + circuit.ccvs.size() * 2 + circuit.voltageSources.size()+circuit.inductors.size(); // 从 DC Voltage Source 之后开始
// for (const auto &vcr : circuit.vcrs)
// {
// int cp = vcr.controlNodePos - 1;
// int cn = vcr.controlNodeNeg - 1;
// int op = vcr.outputNodePos - 1;
// int on = vcr.outputNodeNeg - 1;
// if (on>=0)
// {
// G_temp[on][vcrI_tempndex] += 1;
// G_temp[vcrI_tempndex][on] += 1;
// }
// if (op>=0)
// {
// G_temp[op][vcrI_tempndex] -= 1;
// G_temp[vcrI_tempndex][op] -= 1;
// }
// vcrI_tempndex++;
// }
int vcrI_tempndex = n + circuit.vcvs.size() + circuit.cccs.size() + circuit.ccvs.size() * 2 + circuit.voltageSources.size() + circuit.inductors.size(); // 从 DC Voltage Source 之后开始
for (const auto &vcr : circuit.vcrs)
{
int cp = vcr.controlNodePos - 1;
int cn = vcr.controlNodeNeg - 1;
int op = vcr.outputNodePos - 1;
int on = vcr.outputNodeNeg - 1;
if (op >= 0)
{
G_temp[op][vcrI_tempndex] += 1;
G_temp[vcrI_tempndex][op] += 1;
}
if (on >= 0)
{
G_temp[on][vcrI_tempndex] -= 1;
G_temp[vcrI_tempndex][on] -= 1;
}
vcrI_tempndex++;
}
G = G_temp;
}
@ -842,21 +830,131 @@ void test_vcr(double vi)
circuit.addVoltageSource(1, 0, vi);
circuit.addResistor(1, 2, 1000);
// circuit.addVCCS(2, 0, 1, 0, 1000, "VCR", "LINEAR");
circuit.addVCR(2, 0, 1, 0, 1000, "LINEAR");
circuit.addVCR(2, 0, 2, 1, 100, "LINEAR");
Simulator simulator(circuit);
Vector v = simulator.dc_solve();
std::cout << v[2] <<"\t\t"<< v[3] <<"\t\t"<< v[4] << std::endl;
std::cout << v[2] << "\t\t" << v[3] << "\t\t" << v[4] << std::endl;
}
void test_vcr()
{
for (int i = 0; i < 1500;i++)
{
test_vcr(i*0.001);
}
for (int i = 2; i < 150;i++)
{
test_vcr(i);
}
// test_vcr(10);
// test_vcr(1);
// for (int i = 0; i < 1500; i++)
// {
// test_vcr(i*0.001);
// }
// for (int i = 1; i < 100; i++)
// {
// test_vcr(i);
// }
test_vcr(10);
}
Vector Simulator::dc_solve_new()
{
// OUTPUT: iter_nums, err, V, I
// 注意，前两个元素是迭代次数和误差
Vector V(n + m); // 初始电压矢量
// 牛顿-拉夫森迭代法求解非线性电路
int iter_nums = 0;
double err = 1;
while (err > 1e-8)
{
Matrix G_iter = G; // 复制导纳矩阵
Vector I_iter = I; // 复制电流矢量
// 处理vcr器件
int vcrI_tempndex = n + circuit.vcvs.size() + circuit.cccs.size() + circuit.ccvs.size() * 2 + circuit.voltageSources.size() + circuit.inductors.size(); // 从电感之后开始
for (const auto &vcr : circuit.vcrs)
{
int cp = vcr.controlNodePos - 1;
int cn = vcr.controlNodeNeg - 1;
int op = vcr.outputNodePos - 1;
int on = vcr.outputNodeNeg - 1;
double control_U = (cp >= 0 ? V[cp] : 0) - (cn >= 0 ? V[cn] : 0);
// double out_U = -(op >= 0 ? V[op] : 0) - (on >= 0 ? V[on] : 0);
if (cp >= 0)
{
G_iter[cp][vcrI_tempndex] -= vcr.transresistance * V[vcrI_tempndex];
}
if (cn >= 0)
{
G_iter[cn][vcrI_tempndex] += vcr.transresistance * V[vcrI_tempndex];
}
G_iter[vcrI_tempndex][vcrI_tempndex] -= vcr.transresistance * control_U;
I_iter[vcrI_tempndex] += V[vcrI_tempndex] * vcr.transresistance * control_U;
// 下面的方式使用电阻的方式stamp，但是效果不是很理想
// if (op >= 0)
// {
// G_iter[op][op] += G;
// }
// if (on >= 0)
// {
// G_iter[on][on] += G;
// }
// if (op > 0 && on > 0)
// {
// G_iter[op][on] -= G;
// G_iter[on][op] -= G;
// }
// I_iter[vcrI_tempndex] = out_U;
// vcrI_tempndex++;
}
// 处理二极管
for (const auto &diode : circuit.diodes)
{
int p = diode.nodeP - 1;
int n = diode.nodeN - 1;
double Vd = (p >= 0 ? V[p] : 0) - (n >= 0 ? V[n] : 0);
double Id = diode.Is / (diode.n * 0.02585) * exp(Vd / (diode.n * 0.02585)) * Vd - diode.Is * (exp(Vd / (diode.n * 0.02585)) - 1);
double conductance = diode.Is * exp(Vd / (diode.n * 0.02585)) / (diode.n * 0.02585);
if (p >= 0)
{
I_iter[p] += Id;
G_iter[p][p] += conductance;
}
if (n >= 0)
{
I_iter[n] -= Id;
G_iter[n][n] += conductance;
}
if (p >= 0 && n >= 0)
{
G_iter[p][n] -= conductance;
G_iter[n][p] -= conductance;
}
}
// 解方程
Vector err_V = linear_solve(G_iter, I_iter);
// 检查收敛性
err = (err_V - V).norm();
iter_nums++;
if (iter_nums > 1000)
{
// std::cout << "DC solve failed!" << std::endl;
// std::cout << "Iterations: " << iter_nums << " Error: " << err << std::endl;
break;
}
V = err_V;
}
return V;
}
void test_dio_dcsolver_new(double v)
{
Circuit circuit(3);
circuit.addDiode(1, 0, 1e-12, 1);
circuit.addDiode(0, 1, 1e-12, 1);
circuit.addVoltageSource(2, 0, v);
circuit.addResistor(2, 1, 1000);
Simulator simulator(circuit);
Vector V = simulator.dc_solve_new();
std::cout << v << " " << V[0] << " " << V[1] << " " << V[2] << std::endl;
}
void test_dcsolver_new()
{
for (int i = 1; i < 2000; i++)
{
test_dio_dcsolver_new(sin(i * 0.005));
}
}

6
src/devices.cpp
View File

@ -6,7 +6,7 @@ double get_vcr_G(VCR vcr, double U)
{
if(vcr.transresistance == 0)
{
return 1e20;
return 1e5;
}
return 1/vcr.transresistance;
}
@ -15,11 +15,11 @@ double get_vcr_G(VCR vcr, double U)
double R = vcr.transresistance * U;
if(R == 0)
{
return 1e20;
return 1e5;
}
return 1/R;
}
return 1e20;
return 1e5;
}
double get_voltage(VoltageSource voltage, double time)