245 Reputation

5 years, 137 days

Thank you for the clarification....

@acer Sorry, I missed that point. Thank you for bringing up the difference.

Thank you all for your attention and the ways to make it work!

For the sake of reference, here is the collection of all the answers above.

 > restart;
 > with(LinearAlgebra):
 > with(DynamicSystems):
 > interface(imaginaryunit=j):
 > eq_K1_m4 := K__1 = E__q0*(R__T*E__B*sin(delta) + X__Td*E__B*cos(delta))/(L__aqs*L__l + L__aqs*X__E + L__aqs*L__ads_p + L__l^2 + 2*L__l*X__E + L__l*L__ads_p + R__E^2 + 2*R__E*R__a + R__a^2 + X__E^2 + X__E*L__ads_p) + (X__q - X__dp)*i__q0*(X__Tq*E__B*sin(delta) - R__T*E__B*cos(delta))/(L__aqs*L__l + L__aqs*X__E + L__aqs*L__ads_p + L__l^2 + 2*L__l*X__E + L__l*L__ads_p + R__E^2 + 2*R__E*R__a + R__a^2 + X__E^2 + X__E*L__ads_p);
 (1)
 > eq_K1_m4_desired := K__1 = E__q0*(R__T*E__B*sin(delta) + X__Td*E__B*cos(delta))/Dx + (X__q - X__dp)*i__q0*(X__Tq*E__B*sin(delta) - R__T*E__B*cos(delta))/Dx;
 (2)
 > eq_Dx := Dx = L__aqs*L__l + L__aqs*X__E + L__aqs*L__ads_p + L__l^2 + 2*L__l*X__E + L__l*L__ads_p + R__E^2 + 2*R__E*R__a + R__a^2 + X__E^2 + X__E*L__ads_p;
 (3)
 > denom(op(1, rhs(eq_K1_m4))) - rhs(eq_Dx); # checking to see if the denominator expression is the same as the expression of Dx
 (4)
 > denom(op(2, rhs(eq_K1_m4))) - rhs(eq_Dx); # checking to see if the denominator expression is the same as the expression of Dx
 (5)
 > # 1
 > map2(applyrule, eq_Dx, eq_K1_m4); # did not work.
 (6)
 > # 2
 > subs(eq_Dx, eq_K1_m4); # did not work.
 (7)
 > # 3
 > simplify(eq_K1_m4, {Dx = L__aqs*L__l + L__aqs*X__E + L__aqs*L__ads_p + L__l^2 + 2*L__l*X__E + L__l*L__ads_p + R__E^2 + 2*R__E*R__a + R__a^2 + X__E^2 + X__E*L__ads_p}, [Dx]); # did not work.
 (8)
 > # 4
 > algsubs(eq_Dx, eq_K1_m4); # did not work.
 (9)
 > # 5
 > applyrule(L__aqs*L__l + L__aqs*X__E + L__aqs*L__ads_p + L__l^2 + 2*L__l*X__E + L__l*L__ads_p + R__E^2 + 2*R__E*R__a + R__a^2 + X__E^2 + X__E*L__ads_p = Dx, eq_K1_m4); # did not work.
 (10)
 > Physics:-Substitute((rhs=lhs)(eq_Dx), eq_K1_m4);
 (11)
 > # answer by Carl Love
 > subs(1/rhs(eq_Dx)= 1/Dx, eq_K1_m4);
 (12)
 > subs(denom(rhs(eq_K1_m4))=Dx,eq_K1_m4);
 (13)
 >

For a referece, the following worksheet shows all the answers above:

 > restart;
 > with(LinearAlgebra):
 > with(DynamicSystems):
 > interface(imaginaryunit=j):
 > eq_m1 := m1 = (X__Tq*E__B*sin(delta) - R__T*E__B*cos(delta))/(R__E^2 + 2*R__a*R__E + R__a^2 + X__E^2 + X__E*L__ads*L__fd/(L__fd + L__ads) + 2*X__E*L__l + L__aqs*X__E + L__aqs*L__ads*L__fd/(L__fd + L__ads) + L__aqs*L__l + L__l*L__ads*L__fd/(L__fd + L__ads) + L__l^2);
 (1)
 > eq_m1_desired := m1 = (X__Tq*E__B*sin(delta) - R__T*E__B*cos(delta))/(R__E^2 + 2*R__a*R__E + R__a^2 + X__E^2 + X__E*L__ads_p + 2*X__E*L__l + L__aqs*X__E + L__aqs*L__ads_p + L__aqs*L__l + L__l*L__ads_p + L__l^2);
 (2)
 > ######
 > # by sursumCorda
 (3)
 > # by tomleslie
 (4)
 >

Hello all,

Thank you for all your answers. The attached worksheet is the recap of them.

@Thomas Richard , yes there are different ways of addressing the issue. The last portion of the attached worksheet shows one possbiel way

 > restart;
 > with(LinearAlgebra):
 > interface(imaginaryunit=j):
 > Amat := Matrix(2, 2, [[-0.1428571428*K__D, -0.1081971238], [376.9911185, 0]]);
 (1)
 > A_eig := Eigenvalues(Amat);
 (2)
 > Desired := sqrt((2.000000000*10^(-10))^2 * (1.275510203*10^17*K__D^2 - 1.019733868*10^21));
 (3)
 > # by C_R
 > (identify@expand@factor)~(A_eig);
 (4)
 > # by Acer
 > evalindets(A_eig, &*(float,`+`^(1/2)),u->signum(op(1,u))*(op([2,1],u)*op(1,u)^2)^(1/2));
 (5)
 > # for Thomas Richard
 > restart;
 > with(LinearAlgebra):
 > with(DynamicSystems):
 > interface(imaginaryunit=j):
 > A11 := -K__D / (2*H_gen);
 (6)
 > A12 := -K__S / (2*H_gen);
 (7)
 > A21 := omega__0;
 (8)
 > A22 := 0;
 (9)
 > Amat := <,>;
 (10)
 > Bmat := <1/(2*H_gen), 0>;
 (11)
 > sys_gen := StateSpace(Amat, Bmat);
 > Amat_eigenVals := Eigenvalues(Amat);
 (12)
 > char_poly_gen := CharacteristicPolynomial(sys_gen);
 (13)
 > # Given system parameters
 > Freq := 60;
 (14)
 > sys_MVA_base := MVA * 4;
 (15)
 > H_genx := 3.5;
 (16)
 > omega__0x := evalf(2*Pi*Freq);
 (17)
 > P_gen := 0.9; # in pu
 (18)
 > Q_gen := 0.3; # in pu
 (19)
 > E__t := polar(1.0, convert(36*degrees, radians));
 (20)
 > E__b := polar(0.995, convert(0*degrees, radians));
 (21)
 > Xd__p := j * 0.3;
 (22)
 > Xtr := j * 0.15;
 (23)
 > Xckt1 := j * 0.5;
 (24)
 > Hgen := 3.5;
 (25)
 > ######
 > I__t := conjugate(P_gen +j*Q_gen) / abs(E__t);
 (26)
 > Xd__p * I__t;
 (27)
 > abs(E__t) + Xd__p * I__t;
 (28)
 > E__p := convert(abs(E__t) + Xd__p * I__t, polar);
 (29)
 > temp := convert(argument(E__p), degrees);
 (30)
 > delta__0 := 36*degrees + temp;
 (31)
 > X__T := Xd__p + Xtr + Xckt1;
 (32)
 > # from eq. 12.76
 > K__Sx := (abs(E__p) * abs(E__b) / abs(X__T)) * cos (convert(delta__0, radians));
 (33)
 > Amat_eigenValsNumeric := subs({H_gen=H_genx, K__S=K__Sx, omega__0 = omega__0x}, Amat_eigenVals);
 (34)
 >

I tried.

Thank you!...

 > restart;
 > with(LinearAlgebra):
 > interface(imaginaryunit=j):
 (1)
 (2)
 > eq_5_23 := L__ad__p = evala(rhs(eq_5_23x));
 (3)
 (4)
 > # from Carl Love - 1
 > solve(eq_5_23, {L__ad}); # step 1
 (5)
 > eval(eq_5_22, %); # step 2
 (6)
 > simplify(%); # step 3
 (7)
 > expand(%); # step 4
 (8)
 > expand(simplify(eval(eq_5_22, solve(eq_5_23, {L__ad})))); # single line - combined
 (9)
 > # from Carl Love - 2
 > expand(simplify(eq_5_22, {eq_5_23}, [L__ad])); # using simplify/siderels
 (10)
 > # from Christian Wolinski
 (11)
 >

Sorry!...

@acer You are absolutely right. It does matter!

Thank you!...

@acer Thank you for the answer and the way of using the collect command with no recursion. As noted in the answer to 'tomleslie', it did not matter to me whether the notation 'p' is treated as a function call or a standalone multiplicative factor.

Thank you!...

@tomleslie Thank you for the extra note. The purpose is to check the steps of deriving the final expression, from a series of expressions given in a text book, there is practically no distinction between the two. In other words, if a way becomes preferrable, then it would be utilized for the necessary facilitation.

Yes. the p() with the 'Delta*delta' causes the blockage. Perhaps one way to go around it is to re-write that part from p(Delta*delta) as p*(Delta*delta), which worked with the nested 'collect' command.

Sorry and Thank you...

@Carl Love Thank you for your correction and sorry for my misunderstanding.

Here is the comparison. '1' is your correction, while '2' shows my misunderstanding.

 > ######
 > # 1
 > restart;
 > unprotect(DEtools): DEtools[diffop2de]:= L->     `if`(L::algebraic, `DEtools/diffop2de`, curry(map, procname))(args): protect(DEtools):
 > with(DEtools):
 > _Envdiffopdomain:=[Dt,t]:
 > eq_e5_9 := psi__d0*Delta__delta(t) + psi__q0*Dt*Delta__delta(t)/omega__0 = p_*Delta__psi__d(t)/omega__0 - Delta__psi__q(t);
 (1)
 > diffop2de(eq_e5_9, f(t));
 (2)
 > ######
 > # 2
 > restart;
 > unprotect(DEtools): DEtools[diffop2de]:= L->     `if`(L::algebraic, `DEtools/diffop2de`, curry(map, procname))(args): protect(DEtools):
 > with(DEtools):
 > _Envdiffopdomain:=[Dt,t]:
 > eq_e5_9 := psi__d0*Delta__delta(t) + psi__q0*D*Delta__delta(t)/omega__0 = p_*Delta__psi__d(t)/omega__0 - Delta__psi__q(t);
 (3)
 > diffop2de(eq_e5_9, f(t));
 (4)
 >

Thank you!...

@Carl Love : Thank you for your answer. I was able to re-iterate it by using the attached worksheet (below).

However, still, my question is not fully addressed yet. What I wanted to achieve is the expression labeled as  '(e5.13p) ' in the attached PDF file('Ex5_1_simple.pdf'). Then, if the differentiator operator works as expected, those differentiations associated with constants would go away by themselves.

 > ######
 > # 1
 > restart;
 > unprotect(DEtools): DEtools[diffop2de]:= L->     `if`(L::algebraic, `DEtools/diffop2de`, curry(map, procname))(args): protect(DEtools):
 > with(DEtools):
 > _Envdiffopdomain:=[Dt,t]:
 > eq_e5_9 := psi__d0*Delta__delta(t) + psi__q0*D*Delta__delta(t)/omega__0 = p_*Delta__psi__d(t)/omega__0 - Delta__psi__q(t);
 (1)
 > diffop2de(eq_e5_9, f(t));
 (2)
 > ######
 > # 2
 > restart;
 > unprotect(DEtools): DEtools[diffop2de]:= L->     `if`(L::algebraic, `DEtools/diffop2de`, curry(map, procname))(args): protect(DEtools):
 > with(DEtools):
 > eq_e5_9 := psi__d0*Delta__delta(t) + psi__q0*D*Delta__delta(t)/omega__0 = p_*Delta__psi__d(t)/omega__0 - Delta__psi__q(t);
 (3)
 > diffop2de(eq_e5_9, f(t));
 > ######
 > # 3
 > restart;
 > unprotect(DEtools): DEtools[diffop2de]:= L->     `if`(L::algebraic, `DEtools/diffop2de`, curry(map, procname))(args): protect(DEtools):
 > with(DEtools):
 > _Envdiffopdomain:=[Dx,x]:
 > eq_e5_9 := psi__d0*Delta__delta(t) + psi__q0*D*Delta__delta(t)/omega__0 = p_*Delta__psi__d(t)/omega__0 - Delta__psi__q(t);
 (4)
 > diffop2de(eq_e5_9, f(t));
 (5)
 >

Ex5_1_simple.pdf

Sorry!...

@acer: Thank you for your point. I should have thought that when I posted the question.

Regarding the desired form, please have a look at the following worksheet.

 > ######
 > # E5.13 - 3
 > restart;
 > unprotect(DEtools): DEtools[diffop2de]:= L->     `if`(L::algebraic, `DEtools/diffop2de`, curry(map, procname))(args): protect(DEtools):
 > with(DEtools): _Envdiffopdomain:=[Dt,t]:
 > eq_e5_9 := psi__d0*Delta__delta(t) + psi__q0*D*Delta__delta(t)/omega__0 = D*Delta__psi__d(t)/omega__0 - Delta__psi__q(t);
 (1)
 > diffop2de(eq_e5_9, f(t));
 (2)
 > eq_e5_9_disired := psi__d0*Delta__delta(t) + psi__q0/omega__0*diff(Delta__delta(t),t) = 1/omega__0*diff(Delta__psi__d(t),t) - Delta__psi__q(t);
 (3)
 >

Sorry!...

@acer Thank you for your message and sorry for that missing worksheet. I thought I attached it, but perhaps something went wrong. Here is the worksheet file.

 > ######
 > # E5.13 - 2
 > restart;
 > unprotect(DEtools): DEtools[diffop2de]:= L->     `if`(L::algebraic, `DEtools/diffop2de`, curry(map, procname))(args): protect(DEtools):
 > with(DEtools): _Envdiffopdomain:=[Dt,t]:
 > eq_e5_9 := psi__d0*Delta__delta(t) + psi__q0*D*Delta__delta(t)/omega__0 = D*Delta__psi__d(t)/omega__0 - Delta__psi__q(t);
 (1)
 > diffop2de(eq_e5_9, f(t));
 (2)
 >

Thank you!...

@Carl Love Thank you for the further explanation. Sorry, I should have noticed that the 'unprotect/protect' commands should precede the 'with(DEtools)', Now it works.

 > restart;
 > unprotect(DEtools): DEtools[diffop2de]:= L->     `if`(L::algebraic, `DEtools/diffop2de`, curry(map, procname))(args): protect(DEtools):
 > with(DEtools): _Envdiffopdomain:=[Dt,t]:
 >
 > xxx := psi__d0 + psi__q0*Dt/omega__0 = 1/omega__0 - 1;
 (1)
 > xxx_1 := lhs(xxx) - rhs(xxx);
 (2)
 >
 > xxx_a := diffop2de(xxx, f(t));
 (3)
 > xxx_1_a := diffop2de(xxx_1, f(t));
 (4)
 >

Thank you!...

@Carl Love Thank you for the informative answer. However, I'm afraid the error still persists. The attached worksheet is my attempt with the suggested routine.

 > restart;
 > with(DEtools): _Envdiffopdomain:=[Dt,t]:
 >
 > xxx := psi__d0 + psi__q0*Dt/omega__0 = 1/omega__0 - 1;
 (1)
 > xxx_1 := lhs(xxx) - rhs(xxx);
 (2)
 > unprotect(DEtools): DEtools[diffop2de]:= L->     `if`(L::algebraic, `DEtools/diffop2de`, curry(map, procname))(args): protect(DEtools):
 > xxx_a := diffop2de(xxx, f(t));
 >
 > xxx_1_a := diffop2de(xxx_1, f(t));
 (3)
 >