improve documentation of hydro and steam Turbines and Governos

Signed-off-by: martin.moraga <[email protected]>

dpsim-models/include/dpsim-models/Signal/HydroTurbineGovernor.h

@@ -50,22 +50,25 @@ namespace Signal {
             //Setpoint for mechanical Power (pu)
             Real mPref;
 
-            // Delta Omega = Omega_ref-Omega_meas at k
+            // Delta Omega = Omega_ref-Omega_meas at k and at k-1
             Real mDelOm;
+            Real mDelOm_prev;
 
-            // State Variable of T1 PT1 at k
+            // State Variable of T1 PT1 at k and at k-1
             Real mX1;
-            // State Variable of T1 PT1 at k+1
-            Real mX1_next;
+            Real mX1_prev;
 
-            // State Variable of T3 PT1 at k
+            // State Variable of T3 PT1 at k and at k-1
             Real mX2;
-            // State Variable of T3 PT1 at k+1
-            Real mX2_next;
+            Real mX2_prev;
 
             // The outpur of the Governor at k
             Real mPgv;
 
+            /// Auxiliar variables
+            Real mCa;
+            Real mCb;
+
         public:
             ///
             explicit HydroTurbineGovernor(const String & name) : SimSignalComp(name, name) { }

dpsim-models/include/dpsim-models/Signal/PSS1A.h

@@ -69,6 +69,10 @@ namespace Signal {
 		Real mVs_prev;
 		/// 
 		Real mOmega_prev;
+		///
+		Real mActivePower_prev;
+		///
+		Real mVh_prev;
 
 	protected:
 		/// State Variables
@@ -78,7 +82,7 @@ namespace Signal {
 		Real mV2;
 		/// Output of the second phase compensation block
 		Real mV3;
-		/// PSS output
+		/// PSS output at t=k
 		Real mVs;
 
 	public:

dpsim-models/include/dpsim-models/Signal/SteamTurbineGovernor.h

@@ -21,26 +21,29 @@ namespace Signal {
 
         public:
             /// Droop, the value 1/K in the controller K(1+sT_2)
-            Real R;
+            Real R = 0;
             /// T_1 related to the differentiator in the controlle K(1+sT_2)/(1+sT_1)
-            Real T1;
+            Real T1 = 0;
             /// T_2 related to the differentiator in the controlle K(1+sT_2)
-            Real T2;
+            Real T2 = 0;
             /// Time constant T_3 of the actuator in the Governor
-            Real T3;
+            Real T3 = 0;
 
             // ### Physical constraints ###
             /// Maximum growth rate
-            Real dPmax;
+            Real dPmax = 0;
             /// Minimum decay rate
-            Real dPmin;
+            Real dPmin = 0;
             /// Maximum mechanical power(pu)
-            Real Pmax;
+            Real Pmax = 0;
             /// Minimum mechanical power (pu)
-            Real Pmin;
+            Real Pmin = 0;
 
             /// Setpoint for omega (pu). It is adviced to choose Om_ref=1
-            Real OmRef;
+            Real OmRef = 0;
+
+            /// Proportional gain of anti-Windup
+			Real Kbc = 0;
     };
 
     /// Steam turbine governor, where Omega_ref, P_ref are constant and T_1=0
@@ -58,30 +61,26 @@ namespace Signal {
             /// Setpoint for mechanical Power (pu)
             Real mPref;
 
-            /// ### Variables at time step k-1 ###
-            Real mDelOm_prev;    
-
-            // ### Variables at time step k ###
-            /// Delta Omega = Omega_ref-Omega_meas at k
+            // ### State Variables ###
+            /// Delta Omega = Omega_ref-Omega_meas at t=k, t=k-1 and t=k-2
             Real mDelOm;
-            /// Variable after the rate limiter and before integrator at k
-            Real mDelPgv;
-            /// The outpur of the Governor at k
-            Real mPgv;
-            /// Windup variable
+            Real mDelOm_prev;
+            Real mDelOm_2prev;
+            /// Windup variable at t=k and t=k-1
             Real mP1;
-            /// Windup variable
+            Real mP1_prev;
+            /// p at time t=k-1
             Real mP;
-            /// Windup variable
-            Real mPlim_in;
-
-            // ### Variables at time step k+1 ###
-            /// The outpur of the PT1 with limiters at k+1 (Governor output)
-            Real mPgv_next;
-            ///
-            Real mP1_next;
-            ///
-            Real mPlim_in_next;
+            /// Derivative of Pgv at time k-1
+            Real mDerPgv;
+            /// The output of the Governor at k before the limiter
+            Real mPgvLim;
+            /// The output of the Governor at k after the limiter
+            Real mPgv;
+
+            /// Auxiliar variables
+            Real mCa;
+            Real mCb;
 
         public:
             ///

dpsim-models/src/Base/Base_ReducedOrderSynchronGenerator.cpp

@@ -477,12 +477,6 @@ void Base::ReducedOrderSynchronGenerator<Complex>::initializeFromNodesAndTermina
 	// set initial interface voltage
 	(**mIntfVoltage)(0,0) = mInitVoltage * mBase_V_RMS;
 
-	// set initial value of current
-	(**mIntfCurrent)(0,0) = mInitCurrent * mBase_I_RMS;
-
-	// set initial interface voltage
-	(**mIntfVoltage)(0,0) = mInitVoltage * mBase_V_RMS;
-
 	SPDLOG_LOGGER_DEBUG(this->mSLog,
 		"\n--- Initialization from power flow  ---"
 		"\nInitial Vd (per unit): {:f}"
@@ -526,7 +520,8 @@ template <>
 void Base::ReducedOrderSynchronGenerator<Complex>::mnaCompPreStep(Real time, Int timeStepCount) {
 	mSimTime = time;
 
-	// update governor variables
+	// Calculate mechanical torque at t=k
+	// input is omega at time t=k, output is mech. torque at time t=k
 	if (mHasTurbine && mHasGovernor) {
 		mMechTorque_prev = **mMechTorque;
 		**mMechTorque = mTurbine->step(mGovernor->step(**mOmMech, mTimeStep), mTimeStep);
@@ -536,22 +531,28 @@ void Base::ReducedOrderSynchronGenerator<Complex>::mnaCompPreStep(Real time, Int
 		**mMechTorque = mGovernor->step(**mOmMech, mTimeStep);
 	}
 
-	// calculate mechanical variables at t=k+1 with forward euler
+	// Calculate electrical torque at t=k
 	**mPower = Complex((**mVdq)(0,0) * (**mIdq)(0,0) + (**mVdq)(1,0) * (**mIdq)(1,0),
 					   (**mVdq)(1,0) * (**mIdq)(0,0) - (**mVdq)(0,0) * (**mIdq)(1,0));
 	**mElecTorque = (**mPower).real();
 	//**mElecTorque = (**mPower).real() / **mOmMech;
-	**mOmMech = **mOmMech + mTimeStep * (1. / (2. * mH) * (mMechTorque_prev - **mElecTorque));
-	**mThetaMech = **mThetaMech + mTimeStep * (**mOmMech * mBase_OmMech);
-	**mDelta = **mDelta + mTimeStep * (**mOmMech - 1.) * mBase_OmMech;
-
-	// update exciter and PSS variables
+
+	// calculate v_pss at time k
+	// inputs are: mOmMech, mElecTorque, mVdq at time t=k
 	if (mHasPSS)
 		mVpss = mPSS->step(**mOmMech, **mElecTorque, (**mVdq)(0,0), (**mVdq)(1,0), mTimeStep);
+
+	// calculate e_fd at time k+1, inputs are at t=k
 	if (mHasExciter) {
 		mEf_prev = **mEf;
 		**mEf = mExciter->step((**mVdq)(0,0), (**mVdq)(1,0), mTimeStep, mVpss);		
 	}
+
+	// calculate mechanical variables at t=k+1 with forward euler
+	**mOmMech = **mOmMech + mTimeStep * (1. / (2. * mH) * (mMechTorque_prev - **mElecTorque));
+	//**mOmMech = **mOmMech + mTimeStep * (1. / (2. * mH) * (mMechTorque_prev - **mElecTorque) / **mOmMech);
+	**mThetaMech = **mThetaMech + mTimeStep * (**mOmMech * mBase_OmMech);
+	**mDelta = **mDelta + mTimeStep * (**mOmMech - 1.) * mBase_OmMech;
 
 	stepInPerUnit();
 	(**mRightVector).setZero();
@@ -562,7 +563,8 @@ template <>
 void Base::ReducedOrderSynchronGenerator<Real>::mnaCompPreStep(Real time, Int timeStepCount) {
 	mSimTime = time;
 
-	// update governor variables
+	// Calculate mechanical torque at t=k
+	// input is omega at time t=k, output is mech. torque at time t=k
 	if (mHasTurbine && mHasGovernor) {
 		mMechTorque_prev = **mMechTorque;
 		**mMechTorque = mTurbine->step(mGovernor->step(**mOmMech, mTimeStep), mTimeStep);
@@ -572,22 +574,28 @@ void Base::ReducedOrderSynchronGenerator<Real>::mnaCompPreStep(Real time, Int ti
 		**mMechTorque = mGovernor->step(**mOmMech, mTimeStep);
 	}
 
-	// calculate mechanical variables at t=k+1 with forward euler
+	// Calculate electrical torque at t=k
 	**mPower = Complex((**mVdq0)(0,0) * (**mIdq0)(0,0) + (**mVdq0)(1,0) * (**mIdq0)(1,0),
 					   (**mVdq0)(1,0) * (**mIdq0)(0,0) - (**mVdq0)(0,0) * (**mIdq0)(1,0));
 	**mElecTorque = (**mPower).real();
 	//**mElecTorque = (**mPower).real() / **mOmMech;
-	**mOmMech = **mOmMech + mTimeStep * (1. / (2. * mH) * (mMechTorque_prev - **mElecTorque));
-	**mThetaMech = **mThetaMech + mTimeStep * (**mOmMech * mBase_OmMech);
-	**mDelta = **mDelta + mTimeStep * (**mOmMech - 1.) * mBase_OmMech;
 
-	// update exciter and PSS variables
+	// calculate v_pss at time k
+	// inputs are: mOmMech, mElecTorque, mVdq at time t=k
 	if (mHasPSS)
 		mVpss = mPSS->step(**mOmMech, **mElecTorque, (**mVdq0)(0,0), (**mVdq0)(1,0), mTimeStep);
+
+	// calculate e_fd at time k+1, inputs are at t=k
 	if (mHasExciter) {
 		mEf_prev = **mEf;
 		**mEf = mExciter->step((**mVdq0)(0,0), (**mVdq0)(1,0), mTimeStep, mVpss);	
 	}
+
+	// calculate mechanical variables at t=k+1 with forward euler
+	**mOmMech = **mOmMech + mTimeStep * (1. / (2. * mH) * (mMechTorque_prev - **mElecTorque));
+	//**mOmMech = **mOmMech + mTimeStep * (1. / (2. * mH) * (mMechTorque_prev - **mElecTorque) / **mOmMech);
+	**mThetaMech = **mThetaMech + mTimeStep * (**mOmMech * mBase_OmMech);
+	**mDelta = **mDelta + mTimeStep * (**mOmMech - 1.) * mBase_OmMech;
 
 	stepInPerUnit();
 	(**mRightVector).setZero();

dpsim-models/src/Signal/ExciterStatic.cpp

@@ -66,7 +66,7 @@ void ExciterStatic::initialize(Real Vh_init, Real Ef_init) {
 	mVin = mVe / (mCa + mCb);
 
 	/// Initial reference value
-	mVref = mVin - mVr;
+	mVref = mVin + mVr;
 
 	/// Initial value for auxilar state variable
 	mXb = mVin;
@@ -77,10 +77,10 @@ void ExciterStatic::initialize(Real Vh_init, Real Ef_init) {
 		"\nExciter type: ExciterStatic"
 		"\nInitially set excitation system initial values:"
 		"\ninit Vh: {:e}"
-		"\ninit Ef: {:e}",
+		"\ninit Ef: {:e}"
 		"\nCalculated set poit and auxilary state variables:"
 		"\nVref : {:e}"
-		"\nXb : {:e}",
+		"\nXb : {:e}"
 		"\nVin = {:e}",
 		mVh, mEfd,
 		mVref, mXb, mVin);
@@ -103,15 +103,13 @@ Real ExciterStatic::step(Real mVd, Real mVq, Real dt, Real Vpss) {
 		mVr = mVh;
 
 	/// Compute Vin at time k
-	mVin = mVref - mVr + Vpss;
+	mVin = mVref - mVr_prev + Vpss;
 
 	/// Compute Xb at time k+1 using euler forward
 	mXb = mXb_prev + (mVin - mXb_prev) * dt / mParameters->Tb;
 
-	// Compute Efd at time k+1 using euler forward
-	mVe = mVin * mCa + mXb * mCb - mParameters->Kbc * (mEfd - mEfdLim);
-
-	//Integrator for T_e time constant !with! Wind-up
+	// Compute Efd at time k using euler forward
+	mVe = mVin * mCa + mXb_prev * mCb - mParameters->Kbc * (mEfd - mEfdLim);
 	mEfd = mEfd + (dt / mParameters->Te) * (mParameters->Ka * mVe - mEfd);
 	if (mEfd > mParameters->MaxEfd)
 		mEfdLim = mParameters->MaxEfd;

dpsim-models/src/Signal/HydroTurbine.cpp

@@ -42,9 +42,12 @@ void HydroTurbine::initialize(Real Pminit) {
 Real HydroTurbine::step(Real Pgv, Real dt){
     // the values pre calculated in the previous step are now actual values
     mX1 = mX1_next;
+
     // Caltulating the value of the only state variable for the next step (integration) using the current
-    mX1_next = 2 * dt/mParameters->Tw * (Pgv-mX1) + mX1;
+    mX1_next = 2 * dt / mParameters->Tw * (Pgv - mX1) + mX1;
+
     // Output is the value Pm(k) = -2*Pgv(k)+X1(k) 
     mPm = 3 * mX1 - (2 * Pgv);
+
     return mPm;
 }

dpsim-models/src/Signal/HydroTurbineGovernor.cpp

@@ -28,6 +28,14 @@ void HydroTurbineGovernor::setParameters(std::shared_ptr<Base::GovernorParameter
 	    std::cout << "Type of parameters class of " << this->name() << " has to be HydroGorvernorParameters!" << std::endl;
 	    throw CPS::TypeException();
     }
+    // FIXME: Ein problem entsteht bei der Zerlegung wenn T1=t3 (gleiche polstellen, es muss eine andere Zerlegung benutzt werden)
+    if (mParameters->T1==mParameters->T3) {
+        SPDLOG_LOGGER_ERROR(mSLog, "T1 can not be equal to T3!!!");
+        throw CPS::Exception();
+    }
+
+    mCa = (mParameters->T1 - mParameters->T2) / (mParameters->T1 - mParameters->T3);
+    mCb = (mParameters->T2 - mParameters->T3) / (mParameters->T1 - mParameters->T3);
 }
 
 void HydroTurbineGovernor::initialize(Real Pref){
@@ -36,10 +44,11 @@ void HydroTurbineGovernor::initialize(Real Pref){
         mPgv=Pref;
         mPref=Pref;
         mDelOm = 0;
+        mDelOm_prev = 0;
         mX1 = 0;
-        mX1_next = 0;
+        mX1_prev = 0;
         mX2 = 0;
-        mX2_next = 0;
+        mX2_prev = 0;
         SPDLOG_LOGGER_INFO(mSLog, 
             "Hydro Governor initial values: \n"
 			"\nPref: {:f}"
@@ -57,24 +66,20 @@ void HydroTurbineGovernor::initialize(Real Pref){
 }
 
 Real HydroTurbineGovernor::step(Real Omega, Real dt){
-    // FIXME: Ein problem entsteht bei der Zerlegung wenn T1=t3 (gleiche polstellen, es muss eine andere Zerlegung benutzt werden)
-    // TODO: move cA and cB to initialize() to avoid recomputation in each simulation step
-    const Real cA = (mParameters->T1 - mParameters->T2) / (mParameters->T1 - mParameters->T3);
-    const Real cB = (mParameters->T2 - mParameters->T3) / (mParameters->T1 - mParameters->T3);
-
-    //Calculation of frequency deviation
-    mDelOm=mParameters->OmRef-Omega;
+    // write the values that were calculated in the previous step
+	mDelOm_prev = mDelOm;
+	mX1_prev = mX1;
+    mX2_prev = mX2;
 
-    // Set the values of state variables calculated in last step as atual values for k
-    mX1=mX1_next;
-    mX2=mX2_next;
+   // Calculate the input of the governor for time step k
+	mDelOm = mParameters->OmRef - Omega;    
 
     // Calculation State variables for k+1 with integrators
-    mX1_next= mX1 + dt / mParameters->T1 * (mDelOm - mX1);
-    mX2_next= mX2 + dt / mParameters->T3 * (mDelOm - mX2);
+    mX1 = mX1_prev + dt / mParameters->T1 * (mDelOm_prev - mX1_prev);
+    mX2 = mX2_prev + dt / mParameters->T3 * (mDelOm_prev - mX2_prev);
 
     // Output of the governor before limiter, values from k are used to caltulate output Pgv(k)
-    mPgv = mPref + 1. / mParameters->R * ( mX1 * cA + mX2 * cB);
+    mPgv = mPref + 1. / mParameters->R * ( mX1 * mCa + mX2 * mCb);
 
     if (mPgv>mParameters->Pmax)
         mPgv = mParameters->Pmax;