Confusion on DblPend's output (and maybe SglPend?)

[samm82]

The output requirement of DblPend states that the outputs of the system are $\alpha_1$ and $\alpha_2$,

but the actual output defined in the code (passed to SystemInformation) is just $\theta$ https://github.com/JacquesCarette/Drasil/blob/1acb8ba1553d16d1c949e3b55049fb2242dfe555/code/drasil-example/dblpendulum/lib/Drasil/DblPendulum/Unitals.hs#L35-L36

which shows up in the output table generated (as part of my work on #3259) as:

I'm assuming that the outputs of $\alpha_1$ and $\alpha_2$ are correct, since they are what the IMs output, but changing outputs to be these values (either just the unitals or the IMs themselves) results in the following error which prevents code from being generated:

dblpendulum.EXE: The following outputs cannot be computed: angularAccelerationIM1, angularAccelerationIM2
Unused definitions are: alpha_x1, alpha_y1, p_x1, p_y1, p_x2, p_y2, force
Known values are: l_1, l_2, m_1, m_2, gravitationalMagnitude, pendDisAngle
CallStack (from HasCallStack):
  error, called at lib\Language\Drasil\CodeSpec.hs:179:20 in drasil-code-0.1.9.0-HcWsGH034SJAsZRmMYbtvM:Language.Drasil.CodeSpec
make: *** [Makefile:270: dblpendulum_gen] Error 1

although the SRS output changes to match what is expected

Is there something I'm missing here? How can we change the outputs of the system to be correct so they can actually be computed?

[cd155]

In short, the theta is only for solving the higher-order ODE. It is a hack to make the ODE solver works. It does not belong to the SRS. More like to belong to a design document, and it is very artificial in this case.

We have a discussion discuss it before, you can see more details in this link. https://github.com/JacquesCarette/Drasil/issues/2986#issuecomment-1152504793

[smiths]

Thank you for the information and the link @cd155. There are also additional complications @samm82. I'm not sure how the mistake got in there, but the output of IM1 and IM2 should really be theta_1 and theta_2. Moreover, theta_1 and theta_2 cannot be solved by one IM alone. This is an example of a coupled system of equations. Our simple model of one IM for one output doesn't hold in all cases.

When someone specifies an ODE they aren't usually interested in the rate of change, but the value of the thing that is changing. The double pendulum is a bit confusing, so instead think about Newton's second law: F = m*a. We could rearrange this for a = F/m. This looks like a theory where we are given F and m and return a. However, we don't really care that much about the acceleration, we care about the position as a function of time. We know that a = dv/dt = dp^2/dt^2, where v is the velocity and p is the position. In this case the IM would say given the unbalanced Force (F) over time and the mass, along with the initial position, solve for position as a function of time. The position isn't in the original equation, but that is what we are really interested in finding.

In the pendulum example alpha is the rate of change of the rate of change of the angle. That is, it is the second derivatives of the rotation angle with respect to time.

The output of IM1 and IM2 should be theta_1 and theta_2. I don't actually remember the relationship between theta and theta_1 and theta_2. When the double pendulum code is run (make run RUNARGS=input.txt) we get the time history of theta. @cd155 what is the connection between theta and theta_1 and theta_2. Why don't we output the two angles?

[cd155]

Currently, the Double Pend case study only outputs theta1 in a text file. Although the ODE solver did calculate theta2, but it didn't collect the theta2.

theta is an artifact. Creating it makes it easy to solve ODE equations in ODE solver. In the generated code, the theta is an array which holds four items.

t1  = idx (sy pendDisAngle) (int 0) -- t1 is theta 1
o1  = idx (sy pendDisAngle) (int 1) -- o1 is omega 1
t2  = idx (sy pendDisAngle) (int 2) -- t2 is theta 2
o2  = idx (sy pendDisAngle) (int 3) -- o2 is omega 2

In PDController, we generate the ODE equation for the ODE solver, but, in double pend, we manually write down the ODE equation. The ODE solver requires to put all dependent variables in an array.

Based on IM1 and IM2, we transfer them to a solvable format for the ODE solver, the theta eventually holds four items. Each time the generated code will output a new array. For example, at time 1, it will generate a solution array with [theta1, omega1, theta2, omega2]. By collecting theta1 and theta2, we get the numerical solution for IM1, IM2

I hope this helps clear up some confusion.

[smiths]

Yes, @cd155 that explanation makes sense. The output, in the output file, being called theta confused me. Calling it theta_1 would be better. :smile: Eventually we'll get Drasil to understand higher order ODEs better so that it knows to not just output theta_1, but all the calculated outputs. We could even add alpha_1 and alpha_2 as outputs, if the user wanted these values. (As discussed above, the acceleration output is rarely of interest.)

[JacquesCarette]

@samm82 does this answer your questions?

[samm82]

I think it answers my question in the sense of "I should leave generating the output requirement of DblPendulum to later" 😅

However, I literally just came across SglPendulum: the current output requirement outputs the length of the rod, which is not only an input, but also doesn't come from the linked IM:

My crack at generating this output requirement from the outputs field of SystemInformation gives a more "reasonable" output, but it is still a value that is given as an input ($\theta_p$):

My naive question (i.e., I haven't looked at the rest of SglPendulum enough to answer this, since I have other tasks and probably someone can answer this quicker than me, although I will take a deeper look later if no one else has time either 😅) is should we remove this value from the inputs, should the output be something different, should this output be described in a way that requires Drasil to know more about vectors/ODEs (as for DblPendulum), or should we do something else? This may very well end up as another "leave generating the output requirement for later" situation.

[cd155]

@smiths Should the "Output Values" be θp (displacement angle of the pendulum)? Currently, the single pendulum did not generate any code, but I assume we would like to print θp in the output file.

[smiths]

Yes @cd155, the output should be the displacement angle as a function of time.

[samm82]

Now I've changed the output for SglPend from its unital to its IM, which had the consequence of displaying the IM's name instead of the quantity's (compare its description in the input and output tables below):

Is this OK, or should we change the name of the IM to match the quantity (or be something else entirely)?

[JacquesCarette]

This is a nice illustration of the difference between what theta_p is and what the IM does. Both phrases seem perfect for what they say, but now we need to decide what we really want in our documentation.

[samm82]

A similar thing was uncovered in Projectile (see this comment). However, I think I misunderstood @smiths's comment and used the quantities' terms for the IMs instead of the IMs'; the results of 05f9eb321626d0b4cc86558a03cc35e475e56159 are: Since this table is an aggregation of the quantities, not (necessarily) their IMs, I think it makes sense use this table to use the terms for each quantity, but retain the original IMs' terms when used elsewhere. Thoughts?

[samm82]

This discussion of what we should be calling our IMs will be continued at #3596; for now, the use of the quantities' terms in the IMs (done in Projectile: 05f9eb3) should be reverted, and instead, the generated output table should explicitly use the term of the quantity itself (since this table is really describing the quantities, only using the IMs as their sources).

[smiths]

Yes, the table should use the term of the quantity itself, not the IM for the description. The source column with the IM is a good idea.

[smiths]

@samm82 is there still confusion on the output for single and double pendulum? That is, has the original purpose of the issue been satisfied so that we can close the issue? I know there is work to do in the future to improve double pendulum so that it outputs both theta_1 and theta_2, but that isn't the purpose of this issue.

If your question is answered we should close this issue, or run the risk of the conversation deviating from the original purpose of the issue (which makes finding the discussion in the future almost impossible.)

[samm82]

Good call - I think I understand what the outputs of SglPend are (at least for now), but now that I better understand the next steps for #3259, I really think the IMs of DblPend need to be changed to actually output theta (even if we can only output theta_1 and not theta_2 yet) - this is likely related to #3600. It seems like that should be a separate issue, which will be a blocker for #3259.

[samm82]

I have an idea for a potential workaround that might allow me to "complete" #3259 with DblPend in its current state; then I'll make an issue for removing that workaround once DblPend is updated.

Is #2986 a sufficient issue for tracking the desired changes to the IMs (reworking their outputs), or should I make a new one? @smiths @cd155

[samm82]

Actually, my idea for the workaround won't work without adding a new "secondary" outputs field to SystemInformation for examples with outputs that don't have references, which seems like a big and messy hack. The other option that doesn't involve reworking the existing IMs, as I mentioned in this comment, would be to create an IM for theta.

1. What would that look like?
2. Where is the best place to track that discussion? Here, #2986, or its own separate issue?

[JacquesCarette]

That hack sure seems like too big a hack. Not sure about the answers to any of the questions though.

[smiths]

@samm82, I wonder if the fix you wish to make is simpler than you think? I don't think our integration of IMs and the generated code is that far along. I wouldn't mind if I'm wrong about this, but I feel like we don't actually use the output field in the IM for the generated code? Even if we do use it, the output field for angular acceleration for double pendulum is wrong. Can we just change the output field to be theta? If this is just a display issue, then this should be easy. If the output field has to match the left hand side of the equation, then we have a deeper problem.

What happens if we change the output field for the angular acceleration IM and change it to the angular rotation? From the SRS point of view, we would want to make the following changes:

• IM:calOfAngularAcceleration1 would be renamed IM:calOfAngle1
• Label would change to "Calculation of angle of first rod"
• Output would change to theta_1
• Output Constraint would be empty (I'm not at all sure about the current constraint that alpha_1 > 0. That would imply that the angle would always get bigger in the same direction. The pendulum is swinging back and forth, so the angular acceleration will be both clockwise and counterclockwise.)
• Add a Notes field that says that theta_1 is calculated by solving the ODE in the Equation together with the initial conditions and using IM:calOfAngle2 (which would be a similarly redefined version of IM:calOfAngularAcceleration2).

It doesn't look like we currently mention the initial values that are required for solving the coupled system of ODEs. (Are the initial values somewhere in the SRS @cd155?) It doesn't have to block the current discussion, but we should add initial values to double pendulum.

[cd155]

Initial values didn't show in the SRS, but I manually hard-coded them in the following file.

https://github.com/JacquesCarette/Drasil/blob/2c21ae138a4c65fb2bcf5c52aebcf6fe603263fb/code/drasil-example/dblpend/lib/Drasil/DblPend/ODEs.hs#L20-L22

-- initial time zero unit in second
(exactDbl 0) 
-- final time 20 unit in second
(exactDbl 20)
-- [initial theta1(3*pi/7), initial omega1, initial theta2(3*pi/4), initial omega2] theta(s) unit in radian
[dbl 1.3463968515384828, exactDbl 0, dbl 2.356194490192345, exactDbl 0]

[smiths]

Thanks for the update @cd155. Eventually, the initial values should be added as inputs in the SRS. We should also read these values from a file, rather than hard-code them into the program.

[cd155]

https://github.com/JacquesCarette/Drasil/blob/master/code/datafiles/dblpend/sampleInput.txt

This is the file for some input values in DblPend. We could move the initial values here.

[samm82]

I was able to implement the changes suggested by @smiths in #3610, but changing the outputs field in its SystemInformation to be these IMs in my genOutReqAddSource branch led to the following error (I haven't pushed this change yet due to this error, but I can if that would be helpful to people):

dblpend.EXE: The following outputs cannot be computed: theta_1, theta_2
Unused definitions are: alpha_x1, alpha_y1, p_x1, p_y1, p_x2, p_y2, force
Known values are: l_1, l_2, m_1, m_2, gravitationalMagnitude, pendDisAngle
CallStack (from HasCallStack):
  error, called at lib\Language\Drasil\CodeSpec.hs:180:20 in drasil-code-0.1.9.0-7iTNJMfuwsT5ZCxfZ6jmZv:Language.Drasil.CodeSpec
make: *** [Makefile:270: dblpend_gen] Error 1

My guess as to why this is happening is either:

1. the way to actually "calculate" the thetas is not clear from a code generation perspective,
2. the connection between pendDisAngle and the thetas is not clear, or
3. some combination of the two.

I haven't looked at the DblPend generated code yet, but I'm wondering if @cd155 has any ideas on how to better "link" these thetas to facilitate code generation, especially in light of how these values are currently output to files from this comment above.

[smiths]

@samm82 does the change in #3610 break the generated code, or are you saying that you made an additional change that broke the code? How integrated is the code generation and the IM? Is the output field in the IMs independent of the output list given in SystemInformation?

[samm82]

The change in #3610 doesn't break the generated code, since the outputs of the system are independent of the IM outputs (at least currently). The change I made for generating the output requirement involves actually "linking" these, so now the outputs of the IMs are relevant, which means we likely need to actually implement a way to get theta as an output from the IMs. Hopefully this is feasible, since we already "hardcoded" a way of obtaining these values; from my understanding, now the hard part is linking that up to the IM itself.

[smiths]

Thanks for the explanation @samm82. The fact that the IMs (actually the requirements) aren't connected to the code generator is a long-term thing to fix, but for now at least, it helps us fix the IMs, without breaking the generated code. :smile: We need a way to tell Drasil that we are calculated theta indirectly, by solving an ODE. I think right now it expects the outputs to explicitly be on the left hand side of an equation.

[samm82]

Looking into the code for DblPend, I think we'd want to make the generated code look like this to capture both $\theta_1$ and $\theta_2$ (does this look right @cd155 and @smiths?):

def func_theta(m_1, m_2, L_2, L_1):
    def f(t, theta):
        return [theta[1], (-9.8 * (2.0 * m_1 + m_2) * math.sin(theta[0]) - m_2 * 9.8 * math.sin(theta[0] - 2.0 * theta[2]) - 2.0 * math.sin(theta[0] - theta[2]) * m_2 * (theta[3] ** 2.0 * L_2 + theta[1] ** 2.0 * L_1 * math.cos(theta[0] - theta[2]))) / (L_1 * (2.0 * m_1 + m_2 - m_2 * math.cos(2.0 * theta[0] - 2.0 * theta[2]))), theta[3], 2.0 * math.sin(theta[0] - theta[2]) * (theta[1] ** 2.0 * L_1 * (m_1 + m_2) + 9.8 * (m_1 + m_2) * math.cos(theta[0]) + theta[3] ** 2.0 * L_2 * m_2 * math.cos(theta[0] - theta[2])) / (L_2 * (2.0 * m_1 + m_2 - m_2 * math.cos(2.0 * theta[0] - 2.0 * theta[2])))]

    r = scipy.integrate.ode(f)
    r.set_integrator("dopri5", atol=1.0e-6, rtol=1.0e-6)
    r.set_initial_value([1.3463968515384828, 0.0, 2.356194490192345, 0.0], 0.0)
    theta1 = [[1.3463968515384828, 0.0, 2.356194490192345, 0.0][0]]
    theta2 = [[1.3463968515384828, 0.0, 2.356194490192345, 0.0][2]]
    while r.successful() and r.t < 20.0:
        r.integrate(r.t + 1.0e-3)
        theta1.append(r.y[0])
        theta2.append(r.y[2])

    theta = [theta1, theta2]

    return theta

The current issue where the IMs aren't "linked" to the outputs is because the current output is defined as such: https://github.com/JacquesCarette/Drasil/blob/2039ec288b4acfdcdced106adf843efde53d1f86/code/drasil-example/dblpend/lib/Drasil/DblPend/Unitals.hs#L185-L190

which is then explicitly used when building the ODE: https://github.com/JacquesCarette/Drasil/blob/2039ec288b4acfdcdced106adf843efde53d1f86/code/drasil-example/dblpend/lib/Drasil/DblPend/ODEs.hs#L43-L46

whereas the current outputs of the IMs are built separately from pendDisAngle: https://github.com/JacquesCarette/Drasil/blob/2039ec288b4acfdcdced106adf843efde53d1f86/code/drasil-example/dblpend/lib/Drasil/DblPend/Unitals.hs#L147-L153

Is there a way for us to build a quantity from another like this so the output of an IM could be idx (sy pendDisAngle) (int 0) for example? @JacquesCarette @balacij

[smiths]

@samm82, yes it looks like the code you provided would work to get both theta_1 and theta_2. It took me a long time to confirm that though because the code is strangely written. I guess this is a byproduct of the generator? The theta in the local function isn't the same as the theta in the main function. It would be really nice to have comments, especially one to say that theta is the vector [theta_1, omega_1, theta_2, omega_2] in the local function. It is also strange to say theta1 = [[1.3463968515384828, 0.0, 2.356194490192345, 0.0][0]] instead of just theta1 = 1.3463968515384828. If we are being picky, I'd also like the time in the output, and I'd rather have the theta values for each time value than a vector for theta_1 followed by a vector for theta_2. That is, I'd like

t[0], theta_1[0], theta_2[0]
t[1], theta_1[1], theta_2[1]
....
t[n-1], theta_1[n-1], theta_2[n-1]

Instead of:

theta_1[0], theta_1[1], ... theta_1[n-1]
theta_2[0], theta_2[1], ... theta_2[n-1]

My guess is that I'm making this beyond the scope you intended @samm82. If your goal is to output theta_2, I'm pretty sure the code you gave should work. :smile:

[JacquesCarette]

You're asking for vectors, at least of fixed finite dimensions, to "work properly" when asking for idx (sy pendDisAngle) (int 0). This is a reasonable ask. Unfortunately, I think it only mostly works.

Note that if you're trying 'abuse' vectors (which is often done!) as a poor man's record system... don't.