734{
735 wxFileName fileName(m_data.atpFile);
736 if (!fileName.IsOk()) {
737 errorMsg = wxString::Format(_("Invalid ATP file path for the electromagnetic element \"%s\"."), m_data.name);
738 return false;
739 }
740
741
742 double fundFreq = m_data.frequency;
743 double timeSim = (1.0 / fundFreq) * static_cast<double>(m_data.cyclesToSS);
744
745
746 wxString atpFolder = m_data.atpPath.GetPath();
747
748
749
750 wxExecuteEnv env;
751 env.cwd = atpFolder;
752
753
754 wxTextFile origFile(m_data.atpFile.GetFullPath());
755
756 fileName.SetPath(m_data.atpWorkFolder);
757 wxString name = fileName.GetFullPath();
758 wxTextFile copyFile(fileName.GetFullPath());
759 if (origFile.Open()) {
760 copyFile.Create();
761 for (size_t i = 0; i < origFile.GetLineCount(); i++) {
762 copyFile.AddLine(origFile.GetLine(i));
763 }
764 origFile.Close();
765
766 wxString stepSizeStr = wxString::Format("%.2E", m_data.stepSize);
767 stepSizeStr.Replace(wxT(","), wxT("."));
768 wxString timeStr = wxString::Format("%.6f", timeSim);
769 timeStr.Replace(wxT(","), wxT("."));
770 SetATPParameter(copyFile, wxT("BEGIN NEW DATA CASE"), 1, 0, 8, stepSizeStr);
771 SetATPParameter(copyFile, wxT("BEGIN NEW DATA CASE"), 1, 8, 8, timeStr);
772 SetATPParameter(copyFile, wxT("BEGIN NEW DATA CASE"), 2, 0, 8, wxT("99999999"));
773 SetATPParameter(copyFile, wxT("BEGIN NEW DATA CASE"), 2, 8, 8, wxString::Format("%8d", m_data.recordFrequency));
774
775
776 AddConnectionToNode(copyFile, m_data.atpNodeName);
777
778 copyFile.Write();
779 }
780 else {
781 errorMsg = wxString::Format(_("Fail to open ATP file of the electromagnetic element \"%s\"."), m_data.name);
782 return false;
783 }
784
785
787
788
789
790
791 wxString cmd = m_data.atpPath.GetFullPath() + wxT(" both ") + fileName.GetFullPath() + wxT(" s -R");
792 long pid = wxExecute(cmd, wxEXEC_ASYNC | wxEXEC_HIDE_CONSOLE | wxEXEC_NODISABLE, nullptr, &env);
793
794
795
796
797 if (!bridge.Connect(5000, 100))
798 {
799 if (!CheckLISFile(fileName, errorMsg))
800 {
801 errorMsg = _("Unable to connect to ATP.\n"
802 "Please make sure you are using an ATP version with PSP support.");
803 }
804
805 return false;
806 }
807
808 bridge.SetProcessId(pid);
809
810
811 auto& data = bridge.GetSharedData();
812 data.vrms = std::abs(m_data.puVoltage) * m_data.baseVoltage / std::sqrt(3);
813 data.phase = std::arg(m_data.puVoltage) * 180.0 / M_PI;
814 data.mode = mode;
815
816 m_data.atpSampleData.clear();
817
818
819 while (bridge.WaitATP())
820 {
821 if (!saveRawData || !saveFFTData) {
822 bridge.ReleaseATP();
823 continue;
824 }
825
826 data = bridge.GetSharedData();
827
828
829
830 if (data.stepCount > 0)
831 {
832 const double dt = data.atpStepsize;
833 const double t0 = data.t - (data.stepCount - 1) * dt;
834
835
836 m_data.atpSampleData.reserve(m_data.atpSampleData.size() + data.stepCount);
837
838 std::vector<double> ia, ib, ic;
839
840 if (m_data.useMedianFilter)
841 {
842 ia.reserve(data.stepCount);
843 ib.reserve(data.stepCount);
844 ic.reserve(data.stepCount);
845
846 for (uint32_t i = 0; i < data.stepCount; ++i)
847 {
848 ia.push_back(data.samples[i].current[0]);
849 ib.push_back(data.samples[i].current[1]);
850 ic.push_back(data.samples[i].current[2]);
851 }
852
853 MedianFilter(ia);
854 MedianFilter(ib);
855 MedianFilter(ic);
856 }
857
858 for (uint32_t i = 0; i < data.stepCount; ++i)
859 {
860 const double t = t0 + i * dt;
861
862 const double currentA = m_data.useMedianFilter ? -ia[i] : -data.samples[i].current[0];
863 const double currentB = m_data.useMedianFilter ? -ib[i] : -data.samples[i].current[1];
864 const double currentC = m_data.useMedianFilter ? -ic[i] : -data.samples[i].current[2];
865
866 const double iaMag = std::hypot(-data.phasor.Id[0], -data.phasor.Iq[0]);
867 const double ibMag = std::hypot(-data.phasor.Id[1], -data.phasor.Iq[1]);
868 const double icMag = std::hypot(-data.phasor.Id[2], -data.phasor.Iq[2]);
869 const double iaAng = std::atan2(-data.phasor.Iq[0], -data.phasor.Id[0]);
870 const double ibAng = std::atan2(-data.phasor.Iq[1], -data.phasor.Id[1]);
871 const double icAng = std::atan2(-data.phasor.Iq[2], -data.phasor.Id[2]);
872
873
874 m_data.atpSampleData.push_back({
875 .t = t,
876 .current = { currentA, currentB, currentC },
877 .mag = {iaMag, ibMag, icMag},
878 .angle = {iaAng, ibAng, icAng}
879 });
880 }
881 }
882
883 bridge.ReleaseATP();
884 }
885
886 if (data.terminate < 0)
887 {
888 CheckLISFile(fileName, errorMsg);
889 return false;
890 }
891 constexpr size_t phase = 0;
892
893
894 m_data.current = std::complex<double>(-data.phasor.Id[phase], -data.phasor.Iq[phase]);
895
896
897 wxDir dir(atpFolder);
898 if (dir.IsOpened()) {
899 wxString file;
900 bool cont = dir.GetFirst(&file, wxT("*.tmp"), wxDIR_FILES);
901 while (cont) {
902 wxRemoveFile(atpFolder + wxT("/") + file);
903 cont = dir.GetNext(&file);
904 }
905 }
906
907 if (saveFFTData)
908 {
909
910 double dataStepSize = m_data.atpSampleData[1].t - m_data.atpSampleData[0].t;
911 bool useRemainder = false;
912 size_t sampleCount = m_data.atpSampleData.size();
913 size_t n = static_cast<size_t>(std::ceil(1.0 / (dataStepSize * fundFreq)));
914
915
916
917 if (n > sampleCount)
918 {
919 n = sampleCount;
920 useRemainder = true;
921 }
922
923 double fs = 1.0 / dataStepSize;
924 double df = fs / static_cast<double>(n);
925
926 if (useRemainder) {
927 double rmder = std::abs(std::remainder(fundFreq, df));
928 size_t minRmderN = n;
929
930 while (rmder > 1e-3)
931 {
932 --n;
933 df = fs / static_cast<double>(n);
934 if (std::abs(std::remainder(fundFreq, df)) < rmder)
935 minRmderN = n;
936
937 rmder = std::abs(std::remainder(fundFreq, df));
938 if (n == 0) {
939 n = minRmderN;
940 df = fs / static_cast<double>(n);
941 break;
942 }
943 }
944 }
945
946 double dtWindow = m_data.atpSampleData.back().t - m_data.atpSampleData[sampleCount - n].t;
947
948 fftw_complex* out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * (n / 2 + 1));
949
950 double* in = (double*)fftw_malloc(sizeof(double) * n);
951
952
953 const size_t first = sampleCount - n;
954 for (size_t i = 0; i < n; ++i) {
955 in[i] = m_data.atpSampleData[first + i].current[phase];
956 }
957
958 fftw_plan p = fftw_plan_dft_r2c_1d(n, in, out, FFTW_ESTIMATE);
959
960 fftw_execute(p);
961 fftw_destroy_plan(p);
962
963 double ampCorrection = 2.0 / static_cast<double>(n);
964
965
966 int fundIndex = -1;
967 double fundMagnitude = 0.0;
968
969 for (size_t i = 0; i < (n / 2 + 1); ++i) {
970 double freq = static_cast<double>(i) * df;
971
972 if (std::remainder(freq, fundFreq)) {
973 fundIndex = static_cast<int>(i);
974
975 fundMagnitude =
976 std::hypot(out[i][0], out[i][1]) * ampCorrection;
977
978 break;
979 }
980 }
981
982 if (fundIndex < 0) {
983 errorMsg = wxString::Format(_("Fail to identify the fundamental frequency of the electromagnetic element \"%s\"."), m_data.name);
984
985 fftw_free(in);
986 fftw_free(out);
987
988 return false;
989 }
990
991
992 m_data.currHarmonics.clear();
993
994 for (size_t i = 0; i < (n / 2 + 1); ++i) {
995 double freq = static_cast<double>(i) * df;
996 int order = static_cast<int>(std::round(freq / fundFreq));
997
998 if (order > m_data.numMaxHarmonics)
999 break;
1000
1001 double magnitude = std::hypot(out[i][0], out[i][1]) * ampCorrection;
1002
1003 if ((magnitude / fundMagnitude) > (m_data.harmonicsThreshold / 100.0)) {
1004 m_data.currHarmonics[order] =
1005 std::complex<double>(out[i][0], out[i][1]) *
1006 ampCorrection / std::sqrt(2.0);
1007 }
1008 }
1009
1010 m_data.inFFTData.clear();
1011 m_data.outFFTData.clear();
1012
1013 for (size_t i = 0; i < n; ++i)
1014 {
1015 m_data.inFFTData.emplace_back(m_data.atpSampleData[first + i].t, in[i]);
1016 }
1017
1018 for (size_t i = 0; i < (n / 2 + 1); ++i)
1019 {
1020 std::complex<double> value(out[i][0] * ampCorrection, out[i][1] * ampCorrection);
1021 double freq = static_cast<double>(i) * df;
1022 m_data.outFFTData.emplace_back(freq, value);
1023 }
1024
1025 fftw_free(in);
1026 fftw_free(out);
1027 }
1028
1029 return true;
1030
1031 errorMsg = "Not implemented yet.\n EMTElement.cpp, line 772ish";
1032 return false;
1033
1034 fileName.SetFullName(fileName.GetName() + wxT(".pl4"));
1035 wxFileInputStream pl4File(fileName.GetFullPath());
1036
1037 if (pl4File.IsOk() && pl4File.GetSize() > 0) {
1038 wxDataInputStream store(pl4File);
1039 store.UseBasicPrecisions();
1040 float* buffer = new float[2];
1041 std::vector<double> timeVec;
1042 std::vector<double> valueVec;
1043 double oldTime = 0.0;
1044 int read16Bytes = 0;
1045 while (!pl4File.Eof()) {
1046 store.ReadFloat(buffer, 2);
1047 if (read16Bytes > 11 && (buffer[0] >= oldTime)) {
1048 oldTime = buffer[0];
1049 timeVec.emplace_back(buffer[0]);
1050 valueVec.emplace_back(buffer[1]);
1051 }
1052 read16Bytes++;
1053 }
1054 delete[] buffer;
1055
1056 if (m_data.useMedianFilter)
1057 MedianFilter(valueVec);
1058
1059
1060 double dataStepSize = timeVec[1] - timeVec[0];
1061 bool useRemainder = false;
1062 size_t n = ceil(1.0 / (dataStepSize * fundFreq));
1063
1064
1065 if (n > valueVec.size()) {
1066 n = valueVec.size();
1067 useRemainder = true;
1068 }
1069
1070 double fs = 1.0 / dataStepSize;
1071 double df = fs / static_cast<double>(n);
1072
1073 if (useRemainder) {
1074 double rmder = abs(remainder(fundFreq, df));
1075 int minRmderN = n;
1076 while (rmder > 1e-3) {
1077 n--;
1078 df = fs / static_cast<double>(n);
1079
1080 if (abs(remainder(fundFreq, df)) < rmder) minRmderN = n;
1081
1082 rmder = abs(remainder(fundFreq, df));
1083
1084 if (n <= 0) {
1085 n = minRmderN;
1086 df = fs / static_cast<double>(n);
1087 break;
1088 }
1089 }
1090 }
1091
1092 double dtWindow = timeVec[timeVec.size() - 1] - timeVec[timeVec.size() - n];
1093
1094 fftw_complex* out;
1095 double* in;
1096 fftw_plan p;
1097 out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * (n / 2 + 1));
1098 in = (double*)fftw_malloc(sizeof(double) * n);
1099
1100 for (size_t i = (valueVec.size() - n); i < valueVec.size(); ++i) {
1101 size_t index = i - (valueVec.size() - n);
1102 in[index] = valueVec[i];
1103 }
1104
1105 p = fftw_plan_dft_r2c_1d(n, in, out, FFTW_ESTIMATE);
1106 fftw_execute(p);
1107 fftw_destroy_plan(p);
1108
1109 double ampCorrection = 2.0 / static_cast<double>(n);
1110
1111
1112 int fundIndex = -1;
1113 double fundMagnitude = 0.0;
1114 for (size_t i = 0; i < (n / 2 + 1); i++) {
1115 double freq = static_cast<double>(i) * df;
1116 if (remainder(freq, fundFreq)) {
1117 fundIndex = i;
1118 fundMagnitude = sqrt(out[i][0] * out[i][0] + out[i][1] * out[i][1]) * ampCorrection;
1119 break;
1120 }
1121 }
1122 if (fundIndex < 0) {
1123 errorMsg = wxString::Format(_("Fail to identify the fundamental frequency of the electromagnetic element \"%s\"."), m_data.name);
1124 return false;
1125 }
1126
1127
1128 m_data.currHarmonics.clear();
1129
1130 for (size_t i = 0; i < (n / 2 + 1); i++) {
1131 double freq = static_cast<double>(i) * df;
1132 int order = static_cast<int>(round(freq / fundFreq));
1133
1134 if (order > m_data.numMaxHarmonics) break;
1135
1136 double magnitude = sqrt(out[i][0] * out[i][0] + out[i][1] * out[i][1]) * ampCorrection;
1137 if ((magnitude / fundMagnitude) > (m_data.harmonicsThreshold / 100.0)) {
1138
1139
1140 m_data.currHarmonics[order] = std::complex<double>(out[i][0], out[i][1]) * ampCorrection / sqrt(2.0);
1141 }
1142 }
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161 fftw_free(in);
1162 fftw_free(out);
1163
1164 }
1165 else {
1166 wxString atpErrorMsg = _("No error found.");
1167
1168 fileName.SetFullName(fileName.GetName() + wxT(".lis"));
1169 wxTextFile lisFile(fileName.GetFullPath());
1170
1171 if (!lisFile.Exists()) {
1172 errorMsg = wxString::Format(_("Fail to run ATP file of the electromagnetic element \"%s\".\nThe ATP program did not return any error messages."), m_data.name);
1173 return false;
1174 }
1175 if (lisFile.Open()) {
1176 wxString line = lisFile.GetFirstLine() + "\n";
1177 bool foundError = false;
1178 int lineCount = -1;
1179 atpErrorMsg.Clear();
1180 while (!lisFile.Eof())
1181 {
1182 if (line.Find(wxT("KILL ")) != wxNOT_FOUND) {
1183 foundError = true;
1184 lineCount++;
1185 }
1186 if (line.Find(wxT("----------")) != wxNOT_FOUND) {
1187 foundError = false;
1188 }
1189 if (foundError && lineCount == 1)
1190 atpErrorMsg += line + " ";
1191 line = lisFile.GetNextLine();
1192 }
1193 lisFile.Close();
1194 }
1195 else
1196 {
1197 errorMsg = wxString::Format(_("Fail to run ATP file of the electromagnetic element \"%s\".\nThe ATP program did not return any error messages."), m_data.name);
1198 return false;
1199 }
1200
1201 errorMsg = wxString::Format(_("Fail to run ATP file of the electromagnetic element \"%s\".\nThe ATP returned the following error message:\n\"%s\""), m_data.name, atpErrorMsg);
1202 return false;
1203 }
1204 return true;
1205}
Public Member Functions inherited from Shunt
1.9.8