Python求解多机系统暂态分析

求解多机系统暂态分析

好吧，我不装了，这是一个老师布置的题目，题目也是电力系统分析的C17课后第7题.
先放一下原题：

我想对于这类题，编写代码之前的工作量也很大，先是把等值电路求解出来，然后求发电机的功率特性。我写代码的时候完全是按书本上的分段计算法来写的，所以还蛮容易懂的。
先上代码：

from numpy import *
import matplotlib as plt
import os
import math
from matplotlib import pyplot as pltclass Powerflcal:def __init__(self,h,time_out,simulatetime):self.h = hself.n = int(simulatetime/h+1)self.timeout = time_out# h为步长,n为段数，timeout 为故障切除时间self.gen = [1.5,1,3]#gen1,2,3分别为各发电机的初始有功功率self.delta = [38.64149,29.44408,15.52411,9.19741,23.11738,13.91997]#delta角也是各功率角的初始值,用列表储存并迭代,其中序号0表示delta1，1为2，2为3，3为相对角12，4为13，5为23self.simtime = simulatetimeself.k = [1800*self.h*self.h,18000/7*self.h*self.h,1200*self.h*self.h]self.newp=[1.5,1,3]self.deltap = []self.tempdelta=[]self.deltarecord = mat(zeros((self.n+1,6)))self.precord = mat(zeros((self.n+1,3)))self.precord[0,:] = [0,0,0]self.deltarecord[0,:] = self.delta#这里的数据主要是作为一个记录，k是计算值，newp为每个时间段后的发电机功率，deltap为功率变化量,tempdelta为功率角变化量，deltarecord为矩阵，记录功率角def countpower(self,):#故障状态时计算功率角变化，因为存在着不同状态的切换所以要弄不同的条件for i in linspace(0,self.simtime,self.n):if i == 0:self.scpowercount()self.precord[int(i / self.h), :] = self.newpself.tempdelta = [0.5 * self.k[ii] * self.deltap[ii] for ii in range(len(self.k))]for j in range(len(self.k)):#更新发电机的功率角,为通用公式self.delta[j] = self.delta[j] + self.tempdelta[j]self.delta[3] = self.delta[0] - self.delta[1]self.delta[4] = self.delta[0] - self.delta[2]self.delta[5] = self.delta[1] - self.delta[2]self.deltarecord[int(i/self.h)+1,:] = self.deltaelif i < 0.1:#故障状态时的功率角变化self.scpowercount()self.precord[int(i / self.h), :] = self.newpself.tempdelta = [self.k[ii] * self.deltap[ii]+self.tempdelta[ii] for ii in range(len(self.k))]for j in range(len(self.k)):#更新发电机的功率角,为通用公式self.delta[j] = self.delta[j] + self.tempdelta[j]self.delta[3] = self.delta[0] - self.delta[1]self.delta[4] = self.delta[0] - self.delta[2]self.delta[5] = self.delta[1] - self.delta[2]self.deltarecord[int(i/self.h)+1,:] = self.deltaelif i == 0.1:#切除故障的功率角变化self.cutnowpowercount()self.precord[int(i / self.h), :] = self.newpself.tempdelta = [0.5 * self.k[ii] * self.deltap[ii]+self.tempdelta[ii] for ii in range(len(self.k))]for j in range(len(self.k)):# 更新发电机的功率角,为通用公式self.delta[j] = self.delta[j] + self.tempdelta[j]self.delta[3] = self.delta[0] - self.delta[1]self.delta[4] = self.delta[0] - self.delta[2]self.delta[5] = self.delta[1] - self.delta[2]self.deltarecord[int(i /self.h) + 1, :] = self.deltaelif i> 0.1:self.aftercutpowercount()self.precord[int(i / self.h), :] = self.newpself.tempdelta = [self.k[ii] * self.deltap[ii]+self.tempdelta[ii] for ii in range(len(self.k))]for j in range(len(self.k)):# 更新发电机的功率角,为通用公式self.delta[j] = self.delta[j] + self.tempdelta[j]self.delta[3] = self.delta[0] - self.delta[1]self.delta[4] = self.delta[0] - self.delta[2]self.delta[5] = self.delta[1] - self.delta[2]self.deltarecord[int(i / self.h) + 1, :] = self.deltareturn self.deltarecord,self.precorddef scpowercount(self):#故障状态时发电机功率计算，用故障状态时的功率特性计算self.newp[0] = 0.007093 + 0.509312 * math.sin(math.radians(self.delta[3] + 3.149997)) + 0.402458 * math.sin(math.radians(self.delta[4] + 20.21524))self.newp[1] = 0.110463 + 0.509312* math.sin(math.radians(-self.delta[3] + 3.149997)) + 1.588518 * math.sin(math.radians(self.delta[5] + 20.21523))self.newp[2] = 2.727604 + 0.402458 * math.sin(math.radians(-self.delta[4] + 20.21524)) + 1.588518 * math.sin(math.radians(-self.delta[5] + 20.21523))self.deltap = [self.gen[ii] - self.newp[ii] for ii in range(len(self.gen))]def cutnowpowercount(self):#切除故障时的“功率增量”,即前一状态和后一状态的平均self.newp[0] = 0.007093 + 0.509312 * math.sin(math.radians(self.delta[3] + 3.149997)) + 0.402458 * math.sin(math.radians(self.delta[4] + 20.21524))self.newp[1] = 0.110463 + 0.509312 * math.sin(math.radians(-self.delta[3] + 3.149997)) + 1.588518 * math.sin(math.radians(self.delta[5] + 20.21523))self.newp[2] = 2.727604 + 0.402458 * math.sin(math.radians(-self.delta[4] + 20.21524)) + 1.588518 * math.sin(math.radians(-self.delta[5] + 20.21523))self.newp[0] += 0.04674 + 1.307647*math.sin(math.radians(self.delta[3] + 3.93234)) + 1.033298*math.sin(math.radians(self.delta[4] + 20.99761))self.newp[1] += 0.172054 + 1.307647 * math.sin(math.radians(-self.delta[3] + 3.93234)) + 1.982507 * math.sin(math.radians(self.delta[5] + 20.99761))self.newp[2] += 2.933027 + 1.033298 * math.sin(math.radians(-self.delta[4] + 20.99761)) + 1.982507 * math.sin(math.radians( -self.delta[5] + 20.99761))self.deltap = [2*self.gen[ii] -self.newp[ii] for ii in range(len(self.gen))]def aftercutpowercount(self):self.newp[0] = 0.04674 + 1.307647*math.sin(math.radians(self.delta[3] + 3.93234)) + 1.033298*math.sin(math.radians(self.delta[4] + 20.99761))self.newp[1] = 0.172054 + 1.307647 * math.sin(math.radians(-self.delta[3] + 3.93234)) + 1.982507 * math.sin(math.radians(self.delta[5] + 20.99761))self.newp[2] = 2.933027 + 1.033298 * math.sin(math.radians(-self.delta[4] + 20.99761)) + 1.982507 * math.sin(math.radians(-self.delta[5] + 20.99761))self.deltap = [self.gen[ii] - self.newp[ii] for ii in range(len(self.gen))]print("下面进行潮流计算，只面向课本第17章第7题")
step = input("请输入步长 Δt：")
step = float(step)
timeout = input("请输入仿真截止时间 t：")
timeout = float(timeout)
work3 = Powerflcal(step,0.1,timeout)
# work3 = Powerflcal(0.05,0.1,0.5)
ax1,ax2 = work3.countpower()
print(ax1)
print(ax2)
delta1 = ax1[0:,0].reshape(-1)#转化成行向量
delta1 = delta1.tolist()
delta1 = delta1[0]#转化成列表delta2 = ax1[0:,1].reshape(-1)#转化成行向量
delta2 = delta2.tolist()
delta2 = delta2[0]#转化成列表delta3 = ax1[0:,2].reshape(-1)#转化成行向量
delta3 = delta3.tolist()
delta3 = delta3[0]#转化成列表delta12 = ax1[0:,3].reshape(-1)#转化成行向量
delta12 = delta12.tolist()
delta12 = delta12[0]#转化成列表delta13 = ax1[0:,4].reshape(-1)#转化成行向量
delta13 = delta13.tolist()
delta13 = delta13[0]#转化成列表delta23 = ax1[0:,5].reshape(-1)#转化成行向量
delta23 = delta23.tolist()
delta23 = delta23[0]#转化成列表plt.figure(1)
t=linspace(0,timeout,int(timeout/step+2))
plt.plot(t,delta1,t,delta2,t,delta3,'s-.',MarkerSize = 3)
plt.xlabel('横坐标t',fontproperties = 'simHei',fontsize =20)
plt.ylabel('纵坐标功率角',fontproperties = 'simHei',fontsize =20)
plt.title('三合一图像步长为'+str(step),fontproperties = 'simHei',fontsize =20)
ax = plt.gca()
ax.legend(('$Δ1$','$Δ2$','$Δ3$'),loc = 'lower right',title = 'legend')plt.figure(2)
plt.plot(t,delta12,t,delta13,t,delta23,'s-.',MarkerSize = 3)
plt.xlabel('横坐标t',fontproperties = 'simHei',fontsize =20)
plt.ylabel('纵坐标相对功率角',fontproperties = 'simHei',fontsize =20)
plt.title('三合一图像步长为'+str(step),fontproperties = 'simHei',fontsize =20)
ay = plt.gca()
ay.legend(('$Δ12$','$Δ13$','$Δ23$'),loc = 'lower right',title = 'legend')
plt.show()
os.system("pause")

这些代码都是我自个脑补的，花了一整个上午，原理也很简单，如果慢慢看代码的话就很容易理解，python的类是真的香，self我都是到处用…对了，关于列表和矩阵的用法跟matlab有些许不一样，矩阵是不能直接扩充（c[4,:]=[1,2,3])这种用法是错的，而且append在叠代也不好用，所以要先确定矩阵的维数。

好了，编写完代码就会出现功率角的变化图：

希望对你有所帮助！

Python求解多机系统暂态分析相关推荐

Python 自动贩卖机系统
随着无人新零售经济的崛起,商场.车站.大厦等各种场所都引入了无人饮品自动售货机,方便人们选购自己想要的饮品.购买者选择想要的饮品,通过投币或扫码的方式支付,支付成功后从出货口取出饮品.要求编写代码,利 ...
造纸机系统的全球与中国市场2022-2028年：技术、参与者、趋势、市场规模及占有率研究报告
报告页数: 150 图表数: 100 报告价格:¥16800 本文研究全球与中国市场造纸机系统的发展现状及未来发展趋势,分别从生产和消费的角度分析造纸机系统的主要生产地区.主要消费地区以及主要的生产商 ...
python 物理学中的应用_利用python求解物理学中的双弹簧质能系统详解
前言本文主要给大家介绍了关于利用python求解物理学中双弹簧质能系统的相关内容,分享出来供大家参考学习,下面话不多说了,来一起看看详细的介绍吧. 物理的模型如下: 在这个系统里有两个物体,它们的质 ...
python代码物理_利用python求解物理学中的双弹簧质能系统详解
前言本文主要给大家介绍了关于利用python求解物理学中双弹簧质能系统的相关内容,分享出来供大家参考学习,下面话不多说了,来一起看看详细的介绍吧. 物理的模型如下: 在这个系统里有两个物体,它们的质 ...
运维堡垒机(跳板机)系统 python
相信各位对堡垒机(跳板机)不陌生,为了保证服务器安全,前面加个堡垒机,所有ssh连接都通过堡垒机来完成,堡垒机也需要有身份认证,授权,访问控制,审计等功能,笔者用Python基本实现了上述功能. A ...
mysql求回购率_用户行为分析——回购率、复购率（SQL、Python求解）
有一个多月没有用Python了,有些生疏o(╥﹏╥)o.通过秦路老师的一道题目,分别使用sql和python求解,顺便复习下python点,重点关注[复购率].[回购率]的解法 ☞秦路老师视频讲解(使 ...
Python学习之道-烤机测试日志Log分析统计
Python学习之道-烤机测试日志Log分析统计问题引出一.环境准备二.实践代码 1.初步实现 2.更新CSV文件写入统计结果 3.运行脚本 4.实现遍历多个Log并汇总结果到Excel 三.遇 ...
python atm银行取款系统_python ATM机案例代码
利用目前学的流程控制写的 ''' ATM机需求: 1.登陆输入账号输入密码每日只有3次登陆密码错误的机会,超过3次禁止登陆 2.查询余额 3.存款 4.取款 5.转帐 6.退出 ''' info ...
【人工智能毕设之基于Python+flask+bilstm的评论情感分析系统-哔哩哔哩】 https://b23.tv/QU56eTl
[人工智能毕设之基于Python+flask+bilstm的评论情感分析系统-哔哩哔哩] https://b23.tv/QU56eTl https://b23.tv/QU56eTl

Python求解多机系统暂态分析

求解多机系统暂态分析

Python求解多机系统暂态分析相关推荐

最新文章

热门文章