diff --git a/.gitignore b/.gitignore index a75d72f84..4846a56c8 100644 --- a/.gitignore +++ b/.gitignore @@ -48,6 +48,7 @@ data/settings_v9f.yaml Manifest-v1.10.toml.bak Manifest-v1.11.toml.bak Manifest-v1.11.toml.backup +data/prob*.bin.default data/prob_dynamic_3_seg.bin.default data/prob_dynamic_1.11_3_seg.bin.default data/prob_dynamic_1.10_3_seg.bin.default diff --git a/.gitignore.license b/.gitignore.license deleted file mode 100644 index fa17de081..000000000 --- a/.gitignore.license +++ /dev/null @@ -1,3 +0,0 @@ -SPDX-FileCopyrightText: 2025 Uwe Fechner - -SPDX-License-Identifier: MIT diff --git a/AUTHORS.md b/AUTHORS.md index cc275d19b..e15262068 100644 --- a/AUTHORS.md +++ b/AUTHORS.md @@ -1,6 +1,5 @@ diff --git a/CITATION.cff.license b/CITATION.cff.license deleted file mode 100644 index 48d71b826..000000000 --- a/CITATION.cff.license +++ /dev/null @@ -1,3 +0,0 @@ -SPDX-FileCopyrightText: 2025 Uwe Fechner - -SPDX-License-Identifier: CC-BY-4.0 diff --git a/LICENSES/MPL-2.0.txt b/LICENSES/MPL-2.0.txt new file mode 100644 index 000000000..ee6256cdb --- /dev/null +++ b/LICENSES/MPL-2.0.txt @@ -0,0 +1,373 @@ +Mozilla Public License Version 2.0 +================================== + +1. Definitions +-------------- + +1.1. "Contributor" + means each individual or legal entity that creates, contributes to + the creation of, or owns Covered Software. + +1.2. "Contributor Version" + means the combination of the Contributions of others (if any) used + by a Contributor and that particular Contributor's Contribution. + +1.3. "Contribution" + means Covered Software of a particular Contributor. + +1.4. "Covered Software" + means Source Code Form to which the initial Contributor has attached + the notice in Exhibit A, the Executable Form of such Source Code + Form, and Modifications of such Source Code Form, in each case + including portions thereof. + +1.5. "Incompatible With Secondary Licenses" + means + + (a) that the initial Contributor has attached the notice described + in Exhibit B to the Covered Software; or + + (b) that the Covered Software was made available under the terms of + version 1.1 or earlier of the License, but not also under the + terms of a Secondary License. + +1.6. "Executable Form" + means any form of the work other than Source Code Form. + +1.7. "Larger Work" + means a work that combines Covered Software with other material, in + a separate file or files, that is not Covered Software. + +1.8. "License" + means this document. + +1.9. "Licensable" + means having the right to grant, to the maximum extent possible, + whether at the time of the initial grant or subsequently, any and + all of the rights conveyed by this License. + +1.10. "Modifications" + means any of the following: + + (a) any file in Source Code Form that results from an addition to, + deletion from, or modification of the contents of Covered + Software; or + + (b) any new file in Source Code Form that contains any Covered + Software. + +1.11. "Patent Claims" of a Contributor + means any patent claim(s), including without limitation, method, + process, and apparatus claims, in any patent Licensable by such + Contributor that would be infringed, but for the grant of the + License, by the making, using, selling, offering for sale, having + made, import, or transfer of either its Contributions or its + Contributor Version. + +1.12. "Secondary License" + means either the GNU General Public License, Version 2.0, the GNU + Lesser General Public License, Version 2.1, the GNU Affero General + Public License, Version 3.0, or any later versions of those + licenses. + +1.13. "Source Code Form" + means the form of the work preferred for making modifications. + +1.14. "You" (or "Your") + means an individual or a legal entity exercising rights under this + License. For legal entities, "You" includes any entity that + controls, is controlled by, or is under common control with You. For + purposes of this definition, "control" means (a) the power, direct + or indirect, to cause the direction or management of such entity, + whether by contract or otherwise, or (b) ownership of more than + fifty percent (50%) of the outstanding shares or beneficial + ownership of such entity. + +2. License Grants and Conditions +-------------------------------- + +2.1. Grants + +Each Contributor hereby grants You a world-wide, royalty-free, +non-exclusive license: + +(a) under intellectual property rights (other than patent or trademark) + Licensable by such Contributor to use, reproduce, make available, + modify, display, perform, distribute, and otherwise exploit its + Contributions, either on an unmodified basis, with Modifications, or + as part of a Larger Work; and + +(b) under Patent Claims of such Contributor to make, use, sell, offer + for sale, have made, import, and otherwise transfer either its + Contributions or its Contributor Version. + +2.2. Effective Date + +The licenses granted in Section 2.1 with respect to any Contribution +become effective for each Contribution on the date the Contributor first +distributes such Contribution. + +2.3. Limitations on Grant Scope + +The licenses granted in this Section 2 are the only rights granted under +this License. No additional rights or licenses will be implied from the +distribution or licensing of Covered Software under this License. +Notwithstanding Section 2.1(b) above, no patent license is granted by a +Contributor: + +(a) for any code that a Contributor has removed from Covered Software; + or + +(b) for infringements caused by: (i) Your and any other third party's + modifications of Covered Software, or (ii) the combination of its + Contributions with other software (except as part of its Contributor + Version); or + +(c) under Patent Claims infringed by Covered Software in the absence of + its Contributions. + +This License does not grant any rights in the trademarks, service marks, +or logos of any Contributor (except as may be necessary to comply with +the notice requirements in Section 3.4). + +2.4. Subsequent Licenses + +No Contributor makes additional grants as a result of Your choice to +distribute the Covered Software under a subsequent version of this +License (see Section 10.2) or under the terms of a Secondary License (if +permitted under the terms of Section 3.3). + +2.5. Representation + +Each Contributor represents that the Contributor believes its +Contributions are its original creation(s) or it has sufficient rights +to grant the rights to its Contributions conveyed by this License. + +2.6. Fair Use + +This License is not intended to limit any rights You have under +applicable copyright doctrines of fair use, fair dealing, or other +equivalents. + +2.7. Conditions + +Sections 3.1, 3.2, 3.3, and 3.4 are conditions of the licenses granted +in Section 2.1. + +3. Responsibilities +------------------- + +3.1. Distribution of Source Form + +All distribution of Covered Software in Source Code Form, including any +Modifications that You create or to which You contribute, must be under +the terms of this License. You must inform recipients that the Source +Code Form of the Covered Software is governed by the terms of this +License, and how they can obtain a copy of this License. You may not +attempt to alter or restrict the recipients' rights in the Source Code +Form. + +3.2. Distribution of Executable Form + +If You distribute Covered Software in Executable Form then: + +(a) such Covered Software must also be made available in Source Code + Form, as described in Section 3.1, and You must inform recipients of + the Executable Form how they can obtain a copy of such Source Code + Form by reasonable means in a timely manner, at a charge no more + than the cost of distribution to the recipient; and + +(b) You may distribute such Executable Form under the terms of this + License, or sublicense it under different terms, provided that the + license for the Executable Form does not attempt to limit or alter + the recipients' rights in the Source Code Form under this License. + +3.3. Distribution of a Larger Work + +You may create and distribute a Larger Work under terms of Your choice, +provided that You also comply with the requirements of this License for +the Covered Software. If the Larger Work is a combination of Covered +Software with a work governed by one or more Secondary Licenses, and the +Covered Software is not Incompatible With Secondary Licenses, this +License permits You to additionally distribute such Covered Software +under the terms of such Secondary License(s), so that the recipient of +the Larger Work may, at their option, further distribute the Covered +Software under the terms of either this License or such Secondary +License(s). + +3.4. Notices + +You may not remove or alter the substance of any license notices +(including copyright notices, patent notices, disclaimers of warranty, +or limitations of liability) contained within the Source Code Form of +the Covered Software, except that You may alter any license notices to +the extent required to remedy known factual inaccuracies. + +3.5. Application of Additional Terms + +You may choose to offer, and to charge a fee for, warranty, support, +indemnity or liability obligations to one or more recipients of Covered +Software. However, You may do so only on Your own behalf, and not on +behalf of any Contributor. You must make it absolutely clear that any +such warranty, support, indemnity, or liability obligation is offered by +You alone, and You hereby agree to indemnify every Contributor for any +liability incurred by such Contributor as a result of warranty, support, +indemnity or liability terms You offer. You may include additional +disclaimers of warranty and limitations of liability specific to any +jurisdiction. + +4. Inability to Comply Due to Statute or Regulation +--------------------------------------------------- + +If it is impossible for You to comply with any of the terms of this +License with respect to some or all of the Covered Software due to +statute, judicial order, or regulation then You must: (a) comply with +the terms of this License to the maximum extent possible; and (b) +describe the limitations and the code they affect. Such description must +be placed in a text file included with all distributions of the Covered +Software under this License. Except to the extent prohibited by statute +or regulation, such description must be sufficiently detailed for a +recipient of ordinary skill to be able to understand it. + +5. Termination +-------------- + +5.1. The rights granted under this License will terminate automatically +if You fail to comply with any of its terms. However, if You become +compliant, then the rights granted under this License from a particular +Contributor are reinstated (a) provisionally, unless and until such +Contributor explicitly and finally terminates Your grants, and (b) on an +ongoing basis, if such Contributor fails to notify You of the +non-compliance by some reasonable means prior to 60 days after You have +come back into compliance. Moreover, Your grants from a particular +Contributor are reinstated on an ongoing basis if such Contributor +notifies You of the non-compliance by some reasonable means, this is the +first time You have received notice of non-compliance with this License +from such Contributor, and You become compliant prior to 30 days after +Your receipt of the notice. + +5.2. If You initiate litigation against any entity by asserting a patent +infringement claim (excluding declaratory judgment actions, +counter-claims, and cross-claims) alleging that a Contributor Version +directly or indirectly infringes any patent, then the rights granted to +You by any and all Contributors for the Covered Software under Section +2.1 of this License shall terminate. + +5.3. In the event of termination under Sections 5.1 or 5.2 above, all +end user license agreements (excluding distributors and resellers) which +have been validly granted by You or Your distributors under this License +prior to termination shall survive termination. + +************************************************************************ +* * +* 6. Disclaimer of Warranty * +* ------------------------- * +* * +* Covered Software is provided under this License on an "as is" * +* basis, without warranty of any kind, either expressed, implied, or * +* statutory, including, without limitation, warranties that the * +* Covered Software is free of defects, merchantable, fit for a * +* particular purpose or non-infringing. The entire risk as to the * +* quality and performance of the Covered Software is with You. * +* Should any Covered Software prove defective in any respect, You * +* (not any Contributor) assume the cost of any necessary servicing, * +* repair, or correction. This disclaimer of warranty constitutes an * +* essential part of this License. No use of any Covered Software is * +* authorized under this License except under this disclaimer. * +* * +************************************************************************ + +************************************************************************ +* * +* 7. Limitation of Liability * +* -------------------------- * +* * +* Under no circumstances and under no legal theory, whether tort * +* (including negligence), contract, or otherwise, shall any * +* Contributor, or anyone who distributes Covered Software as * +* permitted above, be liable to You for any direct, indirect, * +* special, incidental, or consequential damages of any character * +* including, without limitation, damages for lost profits, loss of * +* goodwill, work stoppage, computer failure or malfunction, or any * +* and all other commercial damages or losses, even if such party * +* shall have been informed of the possibility of such damages. This * +* limitation of liability shall not apply to liability for death or * +* personal injury resulting from such party's negligence to the * +* extent applicable law prohibits such limitation. Some * +* jurisdictions do not allow the exclusion or limitation of * +* incidental or consequential damages, so this exclusion and * +* limitation may not apply to You. * +* * +************************************************************************ + +8. Litigation +------------- + +Any litigation relating to this License may be brought only in the +courts of a jurisdiction where the defendant maintains its principal +place of business and such litigation shall be governed by laws of that +jurisdiction, without reference to its conflict-of-law provisions. +Nothing in this Section shall prevent a party's ability to bring +cross-claims or counter-claims. + +9. Miscellaneous +---------------- + +This License represents the complete agreement concerning the subject +matter hereof. If any provision of this License is held to be +unenforceable, such provision shall be reformed only to the extent +necessary to make it enforceable. Any law or regulation which provides +that the language of a contract shall be construed against the drafter +shall not be used to construe this License against a Contributor. + +10. Versions of the License +--------------------------- + +10.1. New Versions + +Mozilla Foundation is the license steward. Except as provided in Section +10.3, no one other than the license steward has the right to modify or +publish new versions of this License. Each version will be given a +distinguishing version number. + +10.2. Effect of New Versions + +You may distribute the Covered Software under the terms of the version +of the License under which You originally received the Covered Software, +or under the terms of any subsequent version published by the license +steward. + +10.3. Modified Versions + +If you create software not governed by this License, and you want to +create a new license for such software, you may create and use a +modified version of this License if you rename the license and remove +any references to the name of the license steward (except to note that +such modified license differs from this License). + +10.4. Distributing Source Code Form that is Incompatible With Secondary +Licenses + +If You choose to distribute Source Code Form that is Incompatible With +Secondary Licenses under the terms of this version of the License, the +notice described in Exhibit B of this License must be attached. + +Exhibit A - Source Code Form License Notice +------------------------------------------- + + This Source Code Form is subject to the terms of the Mozilla Public + License, v. 2.0. If a copy of the MPL was not distributed with this + file, You can obtain one at https://mozilla.org/MPL/2.0/. + +If it is not possible or desirable to put the notice in a particular +file, then You may include the notice in a location (such as a LICENSE +file in a relevant directory) where a recipient would be likely to look +for such a notice. + +You may add additional accurate notices of copyright ownership. + +Exhibit B - "Incompatible With Secondary Licenses" Notice +--------------------------------------------------------- + + This Source Code Form is "Incompatible With Secondary Licenses", as + defined by the Mozilla Public License, v. 2.0. diff --git a/Manifest-v1.10.toml.default b/Manifest-v1.10.toml.default index 6021e74cd..47f29efa0 100644 --- a/Manifest-v1.10.toml.default +++ b/Manifest-v1.10.toml.default @@ -2,7 +2,7 @@ julia_version = "1.10.9" manifest_format = "2.0" -project_hash = "c6fc38d04193f2bbac739c68c7e154eab55bd51a" +project_hash = "5a36713063f56fea2f298f969ff86a63bc583017" [[deps.ADTypes]] git-tree-sha1 = "e2478490447631aedba0823d4d7a80b2cc8cdb32" @@ -196,12 +196,6 @@ git-tree-sha1 = "5a97e67919535d6841172016c9530fd69494e5ec" uuid = "2a0fbf3d-bb9c-48f3-b0a9-814d99fd7ab9" version = "0.2.6" -[[deps.CSTParser]] -deps = ["Tokenize"] -git-tree-sha1 = "0157e592151e39fa570645e2b2debcdfb8a0f112" -uuid = "00ebfdb7-1f24-5e51-bd34-a7502290713f" -version = "3.4.3" - [[deps.CSV]] deps = ["CodecZlib", "Dates", "FilePathsBase", "InlineStrings", "Mmap", "Parsers", "PooledArrays", "PrecompileTools", "SentinelArrays", "Tables", "Unicode", "WeakRefStrings", "WorkerUtilities"] git-tree-sha1 = "deddd8725e5e1cc49ee205a1964256043720a6c3" @@ -272,10 +266,10 @@ uuid = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" version = "1.0.3" [[deps.CommonMark]] -deps = ["Crayons", "PrecompileTools"] -git-tree-sha1 = "5fdf00d1979fd4883b44b754fc3423175c9504b4" +deps = ["PrecompileTools"] +git-tree-sha1 = "351d6f4eaf273b753001b2de4dffb8279b100769" uuid = "a80b9123-70ca-4bc0-993e-6e3bcb318db6" -version = "0.8.16" +version = "0.9.1" [[deps.CommonSolve]] git-tree-sha1 = "0eee5eb66b1cf62cd6ad1b460238e60e4b09400c" @@ -882,10 +876,15 @@ uuid = "ae98c720-c025-4a4a-838c-29b094483192" version = "0.2.1" [[deps.JuliaFormatter]] -deps = ["CSTParser", "CommonMark", "DataStructures", "Glob", "PrecompileTools", "TOML", "Tokenize"] -git-tree-sha1 = "59cf7ad64f1b0708a4fa4369879d33bad3239b56" +deps = ["CommonMark", "Glob", "JuliaSyntax", "PrecompileTools", "TOML"] +git-tree-sha1 = "56b382cd34b1a80f63211a0b009461915915bf9e" uuid = "98e50ef6-434e-11e9-1051-2b60c6c9e899" -version = "1.0.62" +version = "2.1.2" + +[[deps.JuliaSyntax]] +git-tree-sha1 = "937da4713526b96ac9a178e2035019d3b78ead4a" +uuid = "70703baa-626e-46a2-a12c-08ffd08c73b4" +version = "0.4.10" [[deps.JumpProcesses]] deps = ["ArrayInterface", "DataStructures", "DiffEqBase", "DiffEqCallbacks", "DocStringExtensions", "FunctionWrappers", "Graphs", "LinearAlgebra", "Markdown", "PoissonRandom", "Random", "RandomNumbers", "RecursiveArrayTools", "Reexport", "SciMLBase", "Setfield", "StaticArrays", "SymbolicIndexingInterface", "UnPack"] @@ -1156,9 +1155,9 @@ version = "0.8.1" [[deps.ModelingToolkit]] deps = ["ADTypes", "AbstractTrees", "ArrayInterface", "BlockArrays", "Combinatorics", "CommonSolve", "Compat", "ConstructionBase", "DataStructures", "DiffEqBase", "DiffEqCallbacks", "DiffEqNoiseProcess", "DiffRules", "DifferentiationInterface", "Distributed", "Distributions", "DocStringExtensions", "DomainSets", "DynamicQuantities", "EnumX", "ExprTools", "FindFirstFunctions", "ForwardDiff", "FunctionWrappers", "FunctionWrappersWrappers", "Graphs", "InteractiveUtils", "JuliaFormatter", "JumpProcesses", "Latexify", "Libdl", "LinearAlgebra", "MLStyle", "Moshi", "NaNMath", "NonlinearSolve", "OffsetArrays", "OrderedCollections", "PrecompileTools", "RecursiveArrayTools", "Reexport", "RuntimeGeneratedFunctions", "SCCNonlinearSolve", "SciMLBase", "SciMLStructures", "Serialization", "Setfield", "SimpleNonlinearSolve", "SparseArrays", "SpecialFunctions", "StaticArrays", "SymbolicIndexingInterface", "SymbolicUtils", "Symbolics", "URIs", "UnPack", "Unitful"] -git-tree-sha1 = "378ca36cc938886760296eb61a0a996cf908fed8" +git-tree-sha1 = "c0a7d2a91d8cce967744918d811d74f30af8cdd6" uuid = "961ee093-0014-501f-94e3-6117800e7a78" -version = "9.72.0" +version = "9.78.0" [deps.ModelingToolkit.extensions] MTKBifurcationKitExt = "BifurcationKit" @@ -1227,9 +1226,9 @@ version = "1.2.0" [[deps.NonlinearSolve]] deps = ["ADTypes", "ArrayInterface", "BracketingNonlinearSolve", "CommonSolve", "ConcreteStructs", "DiffEqBase", "DifferentiationInterface", "FastClosures", "FiniteDiff", "ForwardDiff", "LineSearch", "LinearAlgebra", "LinearSolve", "NonlinearSolveBase", "NonlinearSolveFirstOrder", "NonlinearSolveQuasiNewton", "NonlinearSolveSpectralMethods", "PrecompileTools", "Preferences", "Reexport", "SciMLBase", "SimpleNonlinearSolve", "SparseArrays", "SparseMatrixColorings", "StaticArraysCore", "SymbolicIndexingInterface"] -git-tree-sha1 = "eb1007ce104fc6704ac346859bc53976c44edc38" +git-tree-sha1 = "aeb6fb02e63b4d4f90337ed90ce54ceb4c0efe77" uuid = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" -version = "4.5.0" +version = "4.9.0" [deps.NonlinearSolve.extensions] NonlinearSolveFastLevenbergMarquardtExt = "FastLevenbergMarquardt" @@ -1943,9 +1942,9 @@ version = "0.2.2" [[deps.SymbolicUtils]] deps = ["AbstractTrees", "ArrayInterface", "Bijections", "ChainRulesCore", "Combinatorics", "ConstructionBase", "DataStructures", "DocStringExtensions", "DynamicPolynomials", "ExproniconLite", "LinearAlgebra", "MultivariatePolynomials", "NaNMath", "Setfield", "SparseArrays", "SpecialFunctions", "StaticArrays", "SymbolicIndexingInterface", "TaskLocalValues", "TermInterface", "TimerOutputs", "Unityper", "WeakValueDicts"] -git-tree-sha1 = "ca5929df933a8b7272bc7f01dcd14b2d976c56e2" +git-tree-sha1 = "fa63e8f55e99aee528951ba26544403b09645979" uuid = "d1185830-fcd6-423d-90d6-eec64667417b" -version = "3.25.1" +version = "3.29.0" [deps.SymbolicUtils.extensions] SymbolicUtilsLabelledArraysExt = "LabelledArrays" @@ -2064,11 +2063,6 @@ git-tree-sha1 = "fe7046d2b5bc1d31cde8fd19fad7c5506e3960b4" uuid = "21f18d07-b854-4dab-86f0-c15a3821819a" version = "0.1.5" -[[deps.Tokenize]] -git-tree-sha1 = "468b4685af4abe0e9fd4d7bf495a6554a6276e75" -uuid = "0796e94c-ce3b-5d07-9a54-7f471281c624" -version = "0.5.29" - [[deps.TranscodingStreams]] git-tree-sha1 = "0c45878dcfdcfa8480052b6ab162cdd138781742" uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa" @@ -2133,9 +2127,9 @@ version = "0.21.71" [[deps.VortexStepMethod]] deps = ["Colors", "DefaultApplication", "DelimitedFiles", "DifferentiationInterface", "FiniteDiff", "Interpolations", "LaTeXStrings", "LinearAlgebra", "Logging", "Measures", "NonlinearSolve", "Parameters", "Pkg", "PreallocationTools", "PrecompileTools", "SciMLBase", "Serialization", "StaticArrays", "Statistics", "Timers", "Xfoil", "YAML"] -git-tree-sha1 = "fd7616ebd3dc9fc2c0e5e6abf7ba641aa52be444" +git-tree-sha1 = "a0fa008357206434acfa2acc588b84764478df62" uuid = "ed3cd733-9f0f-46a9-93e0-89b8d4998dd9" -version = "1.2.5" +version = "1.2.6" [deps.VortexStepMethod.extensions] VortexStepMethodControlPlotsExt = "ControlPlots" diff --git a/Manifest-v1.10.toml.default.license b/Manifest-v1.10.toml.default.license deleted file mode 100644 index c52bffc3b..000000000 --- a/Manifest-v1.10.toml.default.license +++ /dev/null @@ -1,2 +0,0 @@ -SPDX-FileCopyrightText: 2025 Uwe Fechner -SPDX-License-Identifier: MIT diff --git a/Manifest-v1.10.toml.license b/Manifest-v1.10.toml.license deleted file mode 100644 index fa17de081..000000000 --- a/Manifest-v1.10.toml.license +++ /dev/null @@ -1,3 +0,0 @@ -SPDX-FileCopyrightText: 2025 Uwe Fechner - -SPDX-License-Identifier: MIT diff --git a/Manifest-v1.11.toml.default b/Manifest-v1.11.toml.default index 74d5732c0..d2fd1786c 100644 --- a/Manifest-v1.11.toml.default +++ b/Manifest-v1.11.toml.default @@ -2,7 +2,7 @@ julia_version = "1.11.5" manifest_format = "2.0" -project_hash = "59d5d1e07c5eb8899b6018cabcf39c8e0944f5c2" +project_hash = "f79cb96ee087c4eaf524e3c1be0490737e6246ce" [[deps.ADTypes]] git-tree-sha1 = "e2478490447631aedba0823d4d7a80b2cc8cdb32" @@ -178,12 +178,13 @@ version = "1.6.3" [[deps.BracketingNonlinearSolve]] deps = ["CommonSolve", "ConcreteStructs", "NonlinearSolveBase", "PrecompileTools", "Reexport", "SciMLBase"] -git-tree-sha1 = "637ebe439ba587828fd997b7810d8171eed2ea1b" +git-tree-sha1 = "a9014924595b7a2c1dd14aac516e38fa10ada656" uuid = "70df07ce-3d50-431d-a3e7-ca6ddb60ac1e" -version = "1.2.0" -weakdeps = ["ForwardDiff"] +version = "1.3.0" +weakdeps = ["ChainRulesCore", "ForwardDiff"] [deps.BracketingNonlinearSolve.extensions] + BracketingNonlinearSolveChainRulesCoreExt = ["ChainRulesCore", "ForwardDiff"] BracketingNonlinearSolveForwardDiffExt = "ForwardDiff" [[deps.CEnum]] @@ -197,12 +198,6 @@ git-tree-sha1 = "5a97e67919535d6841172016c9530fd69494e5ec" uuid = "2a0fbf3d-bb9c-48f3-b0a9-814d99fd7ab9" version = "0.2.6" -[[deps.CSTParser]] -deps = ["Tokenize"] -git-tree-sha1 = "0157e592151e39fa570645e2b2debcdfb8a0f112" -uuid = "00ebfdb7-1f24-5e51-bd34-a7502290713f" -version = "3.4.3" - [[deps.CSV]] deps = ["CodecZlib", "Dates", "FilePathsBase", "InlineStrings", "Mmap", "Parsers", "PooledArrays", "PrecompileTools", "SentinelArrays", "Tables", "Unicode", "WeakRefStrings", "WorkerUtilities"] git-tree-sha1 = "deddd8725e5e1cc49ee205a1964256043720a6c3" @@ -271,10 +266,10 @@ uuid = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" version = "1.0.3" [[deps.CommonMark]] -deps = ["Crayons", "PrecompileTools"] -git-tree-sha1 = "5fdf00d1979fd4883b44b754fc3423175c9504b4" +deps = ["PrecompileTools"] +git-tree-sha1 = "351d6f4eaf273b753001b2de4dffb8279b100769" uuid = "a80b9123-70ca-4bc0-993e-6e3bcb318db6" -version = "0.8.16" +version = "0.9.1" [[deps.CommonSolve]] git-tree-sha1 = "0eee5eb66b1cf62cd6ad1b460238e60e4b09400c" @@ -618,9 +613,9 @@ uuid = "4e289a0a-7415-4d19-859d-a7e5c4648b56" version = "1.0.5" [[deps.EnzymeCore]] -git-tree-sha1 = "1eb59f40a772d0fbd4cb75e00b3fa7f5f79c975a" +git-tree-sha1 = "7d7822a643c33bbff4eab9c87ca8459d7c688db0" uuid = "f151be2c-9106-41f4-ab19-57ee4f262869" -version = "0.8.9" +version = "0.8.11" weakdeps = ["Adapt"] [deps.EnzymeCore.extensions] @@ -886,16 +881,21 @@ uuid = "ae98c720-c025-4a4a-838c-29b094483192" version = "0.2.1" [[deps.JuliaFormatter]] -deps = ["CSTParser", "CommonMark", "DataStructures", "Glob", "PrecompileTools", "TOML", "Tokenize"] -git-tree-sha1 = "59cf7ad64f1b0708a4fa4369879d33bad3239b56" +deps = ["CommonMark", "Glob", "JuliaSyntax", "PrecompileTools", "TOML"] +git-tree-sha1 = "56b382cd34b1a80f63211a0b009461915915bf9e" uuid = "98e50ef6-434e-11e9-1051-2b60c6c9e899" -version = "1.0.62" +version = "2.1.2" + +[[deps.JuliaSyntax]] +git-tree-sha1 = "937da4713526b96ac9a178e2035019d3b78ead4a" +uuid = "70703baa-626e-46a2-a12c-08ffd08c73b4" +version = "0.4.10" [[deps.JumpProcesses]] -deps = ["ArrayInterface", "DataStructures", "DiffEqBase", "DocStringExtensions", "FunctionWrappers", "Graphs", "LinearAlgebra", "Markdown", "PoissonRandom", "Random", "RandomNumbers", "RecursiveArrayTools", "Reexport", "SciMLBase", "Setfield", "StaticArrays", "SymbolicIndexingInterface", "UnPack"] -git-tree-sha1 = "f2bdec5b4580414aee3178c8caa6e46c344c0bbc" +deps = ["ArrayInterface", "DataStructures", "DiffEqBase", "DiffEqCallbacks", "DocStringExtensions", "FunctionWrappers", "Graphs", "LinearAlgebra", "Markdown", "PoissonRandom", "Random", "RandomNumbers", "RecursiveArrayTools", "Reexport", "SciMLBase", "Setfield", "StaticArrays", "SymbolicIndexingInterface", "UnPack"] +git-tree-sha1 = "216c196df09c8b80a40a2befcb95760eb979bcfd" uuid = "ccbc3e58-028d-4f4c-8cd5-9ae44345cda5" -version = "9.14.3" +version = "9.15.0" weakdeps = ["FastBroadcast"] [[deps.KLU]] @@ -929,19 +929,21 @@ version = "1.4.0" [[deps.Latexify]] deps = ["Format", "InteractiveUtils", "LaTeXStrings", "MacroTools", "Markdown", "OrderedCollections", "Requires"] -git-tree-sha1 = "cd10d2cc78d34c0e2a3a36420ab607b611debfbb" +git-tree-sha1 = "4f34eaabe49ecb3fb0d58d6015e32fd31a733199" uuid = "23fbe1c1-3f47-55db-b15f-69d7ec21a316" -version = "0.16.7" +version = "0.16.8" [deps.Latexify.extensions] DataFramesExt = "DataFrames" SparseArraysExt = "SparseArrays" SymEngineExt = "SymEngine" + TectonicExt = "tectonic_jll" [deps.Latexify.weakdeps] DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" SymEngine = "123dc426-2d89-5057-bbad-38513e3affd8" + tectonic_jll = "d7dd28d6-a5e6-559c-9131-7eb760cdacc5" [[deps.LayoutPointers]] deps = ["ArrayInterface", "LinearAlgebra", "ManualMemory", "SIMDTypes", "Static", "StaticArrayInterface"] @@ -1165,9 +1167,9 @@ version = "0.8.1" [[deps.ModelingToolkit]] deps = ["ADTypes", "AbstractTrees", "ArrayInterface", "BlockArrays", "Combinatorics", "CommonSolve", "Compat", "ConstructionBase", "DataStructures", "DiffEqBase", "DiffEqCallbacks", "DiffEqNoiseProcess", "DiffRules", "DifferentiationInterface", "Distributed", "Distributions", "DocStringExtensions", "DomainSets", "DynamicQuantities", "EnumX", "ExprTools", "FindFirstFunctions", "ForwardDiff", "FunctionWrappers", "FunctionWrappersWrappers", "Graphs", "InteractiveUtils", "JuliaFormatter", "JumpProcesses", "Latexify", "Libdl", "LinearAlgebra", "MLStyle", "Moshi", "NaNMath", "NonlinearSolve", "OffsetArrays", "OrderedCollections", "PrecompileTools", "RecursiveArrayTools", "Reexport", "RuntimeGeneratedFunctions", "SCCNonlinearSolve", "SciMLBase", "SciMLStructures", "Serialization", "Setfield", "SimpleNonlinearSolve", "SparseArrays", "SpecialFunctions", "StaticArrays", "SymbolicIndexingInterface", "SymbolicUtils", "Symbolics", "URIs", "UnPack", "Unitful"] -git-tree-sha1 = "378ca36cc938886760296eb61a0a996cf908fed8" +git-tree-sha1 = "c0a7d2a91d8cce967744918d811d74f30af8cdd6" uuid = "961ee093-0014-501f-94e3-6117800e7a78" -version = "9.72.0" +version = "9.78.0" [deps.ModelingToolkit.extensions] MTKBifurcationKitExt = "BifurcationKit" @@ -1236,9 +1238,9 @@ version = "1.2.0" [[deps.NonlinearSolve]] deps = ["ADTypes", "ArrayInterface", "BracketingNonlinearSolve", "CommonSolve", "ConcreteStructs", "DiffEqBase", "DifferentiationInterface", "FastClosures", "FiniteDiff", "ForwardDiff", "LineSearch", "LinearAlgebra", "LinearSolve", "NonlinearSolveBase", "NonlinearSolveFirstOrder", "NonlinearSolveQuasiNewton", "NonlinearSolveSpectralMethods", "PrecompileTools", "Preferences", "Reexport", "SciMLBase", "SimpleNonlinearSolve", "SparseArrays", "SparseMatrixColorings", "StaticArraysCore", "SymbolicIndexingInterface"] -git-tree-sha1 = "eb1007ce104fc6704ac346859bc53976c44edc38" +git-tree-sha1 = "aeb6fb02e63b4d4f90337ed90ce54ceb4c0efe77" uuid = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" -version = "4.5.0" +version = "4.9.0" [deps.NonlinearSolve.extensions] NonlinearSolveFastLevenbergMarquardtExt = "FastLevenbergMarquardt" @@ -1360,9 +1362,9 @@ version = "1.8.1" [[deps.OrdinaryDiffEqBDF]] deps = ["ADTypes", "ArrayInterface", "DiffEqBase", "FastBroadcast", "LinearAlgebra", "MacroTools", "MuladdMacro", "OrdinaryDiffEqCore", "OrdinaryDiffEqDifferentiation", "OrdinaryDiffEqNonlinearSolve", "OrdinaryDiffEqSDIRK", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "StaticArrays", "TruncatedStacktraces"] -git-tree-sha1 = "a72bf554d5fd1f33a8d2aead3562eddd28ba4c76" +git-tree-sha1 = "42755bd13fe56e9d9ce1bc005f8b206a6b56b731" uuid = "6ad6398a-0878-4a85-9266-38940aa047c8" -version = "1.5.0" +version = "1.5.1" [[deps.OrdinaryDiffEqCore]] deps = ["ADTypes", "Accessors", "Adapt", "ArrayInterface", "DataStructures", "DiffEqBase", "DocStringExtensions", "EnumX", "FastBroadcast", "FastClosures", "FastPower", "FillArrays", "FunctionWrappersWrappers", "InteractiveUtils", "LinearAlgebra", "Logging", "MacroTools", "MuladdMacro", "Polyester", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLOperators", "SciMLStructures", "SimpleUnPack", "Static", "StaticArrayInterface", "StaticArraysCore", "SymbolicIndexingInterface", "TruncatedStacktraces"] @@ -1376,15 +1378,15 @@ weakdeps = ["EnzymeCore"] [[deps.OrdinaryDiffEqDifferentiation]] deps = ["ADTypes", "ArrayInterface", "ConcreteStructs", "ConstructionBase", "DiffEqBase", "DifferentiationInterface", "FastBroadcast", "FiniteDiff", "ForwardDiff", "FunctionWrappersWrappers", "LinearAlgebra", "LinearSolve", "OrdinaryDiffEqCore", "SciMLBase", "SciMLOperators", "SparseArrays", "SparseMatrixColorings", "StaticArrayInterface", "StaticArrays"] -git-tree-sha1 = "315d25dd06614e199973cc13d22e533073bd7458" +git-tree-sha1 = "c78060115fa4ea9d70ac47fa49496acbc630aefa" uuid = "4302a76b-040a-498a-8c04-15b101fed76b" -version = "1.9.0" +version = "1.9.1" [[deps.OrdinaryDiffEqNonlinearSolve]] deps = ["ADTypes", "ArrayInterface", "DiffEqBase", "FastBroadcast", "FastClosures", "ForwardDiff", "LinearAlgebra", "LinearSolve", "MuladdMacro", "NonlinearSolve", "OrdinaryDiffEqCore", "OrdinaryDiffEqDifferentiation", "PreallocationTools", "RecursiveArrayTools", "SciMLBase", "SciMLOperators", "SciMLStructures", "SimpleNonlinearSolve", "StaticArrays"] -git-tree-sha1 = "2f956f14c97ff507e855703ac760d513f7c3e372" +git-tree-sha1 = "ffdb0f5207b0e30f8b1edf99b3b9546d9c48ccaf" uuid = "127b3ac7-2247-4354-8eb6-78cf4e7c58e8" -version = "1.9.0" +version = "1.10.0" [[deps.OrdinaryDiffEqSDIRK]] deps = ["ADTypes", "DiffEqBase", "FastBroadcast", "LinearAlgebra", "MacroTools", "MuladdMacro", "OrdinaryDiffEqCore", "OrdinaryDiffEqDifferentiation", "OrdinaryDiffEqNonlinearSolve", "RecursiveArrayTools", "Reexport", "SciMLBase", "TruncatedStacktraces"] @@ -1427,9 +1429,9 @@ version = "0.4.4" [[deps.Polyester]] deps = ["ArrayInterface", "BitTwiddlingConvenienceFunctions", "CPUSummary", "IfElse", "ManualMemory", "PolyesterWeave", "Static", "StaticArrayInterface", "StrideArraysCore", "ThreadingUtilities"] -git-tree-sha1 = "2082cc4be5e765bd982ed04ea06c068f4f702410" +git-tree-sha1 = "6f7cd22a802094d239824c57d94c8e2d0f7cfc7d" uuid = "f517fe37-dbe3-4b94-8317-1923a5111588" -version = "0.7.17" +version = "0.7.18" [[deps.PolyesterWeave]] deps = ["BitTwiddlingConvenienceFunctions", "CPUSummary", "IfElse", "Static", "ThreadingUtilities"] @@ -1640,9 +1642,9 @@ weakdeps = ["RecipesBase"] [[deps.RuntimeGeneratedFunctions]] deps = ["ExprTools", "SHA", "Serialization"] -git-tree-sha1 = "7cb9d10026d630ce2dd2a1fc6006a3d5041b34c0" +git-tree-sha1 = "86a8a8b783481e1ea6b9c91dd949cb32191f8ab4" uuid = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47" -version = "0.5.14" +version = "0.5.15" [[deps.SCCNonlinearSolve]] deps = ["CommonSolve", "PrecompileTools", "Reexport", "SciMLBase", "SymbolicIndexingInterface"] @@ -1861,9 +1863,9 @@ weakdeps = ["SparseArrays"] [[deps.StatsAPI]] deps = ["LinearAlgebra"] -git-tree-sha1 = "1ff449ad350c9c4cbc756624d6f8a8c3ef56d3ed" +git-tree-sha1 = "9d72a13a3f4dd3795a195ac5a44d7d6ff5f552ff" uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0" -version = "1.7.0" +version = "1.7.1" [[deps.StatsBase]] deps = ["AliasTables", "DataAPI", "DataStructures", "LinearAlgebra", "LogExpFunctions", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"] @@ -1971,9 +1973,9 @@ version = "0.2.2" [[deps.SymbolicUtils]] deps = ["AbstractTrees", "ArrayInterface", "Bijections", "ChainRulesCore", "Combinatorics", "ConstructionBase", "DataStructures", "DocStringExtensions", "DynamicPolynomials", "ExproniconLite", "LinearAlgebra", "MultivariatePolynomials", "NaNMath", "Setfield", "SparseArrays", "SpecialFunctions", "StaticArrays", "SymbolicIndexingInterface", "TaskLocalValues", "TermInterface", "TimerOutputs", "Unityper", "WeakValueDicts"] -git-tree-sha1 = "ca5929df933a8b7272bc7f01dcd14b2d976c56e2" +git-tree-sha1 = "fa63e8f55e99aee528951ba26544403b09645979" uuid = "d1185830-fcd6-423d-90d6-eec64667417b" -version = "3.25.1" +version = "3.29.0" [deps.SymbolicUtils.extensions] SymbolicUtilsLabelledArraysExt = "LabelledArrays" @@ -1985,9 +1987,9 @@ version = "3.25.1" [[deps.Symbolics]] deps = ["ADTypes", "ArrayInterface", "Bijections", "CommonWorldInvalidations", "ConstructionBase", "DataStructures", "DiffRules", "Distributions", "DocStringExtensions", "DomainSets", "DynamicPolynomials", "LaTeXStrings", "Latexify", "Libdl", "LinearAlgebra", "LogExpFunctions", "MacroTools", "Markdown", "NaNMath", "OffsetArrays", "PrecompileTools", "Primes", "RecipesBase", "Reexport", "RuntimeGeneratedFunctions", "SciMLBase", "Setfield", "SparseArrays", "SpecialFunctions", "StaticArraysCore", "SymbolicIndexingInterface", "SymbolicLimits", "SymbolicUtils", "TermInterface"] -git-tree-sha1 = "c8bc71edb71a1dcbb241a75d5f57f40804493678" +git-tree-sha1 = "e14834f421edaa8a30493f7864dfc8582855bb3c" uuid = "0c5d862f-8b57-4792-8d23-62f2024744c7" -version = "6.39.1" +version = "6.40.0" [deps.Symbolics.extensions] SymbolicsForwardDiffExt = "ForwardDiff" @@ -2024,9 +2026,9 @@ version = "1.0.1" [[deps.Tables]] deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "OrderedCollections", "TableTraits"] -git-tree-sha1 = "598cd7c1f68d1e205689b1c2fe65a9f85846f297" +git-tree-sha1 = "f2c1efbc8f3a609aadf318094f8fc5204bdaf344" uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" -version = "1.12.0" +version = "1.12.1" [[deps.Tar]] deps = ["ArgTools", "SHA"] @@ -2093,11 +2095,6 @@ git-tree-sha1 = "fe7046d2b5bc1d31cde8fd19fad7c5506e3960b4" uuid = "21f18d07-b854-4dab-86f0-c15a3821819a" version = "0.1.5" -[[deps.Tokenize]] -git-tree-sha1 = "468b4685af4abe0e9fd4d7bf495a6554a6276e75" -uuid = "0796e94c-ce3b-5d07-9a54-7f471281c624" -version = "0.5.29" - [[deps.TranscodingStreams]] git-tree-sha1 = "0c45878dcfdcfa8480052b6ab162cdd138781742" uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa" @@ -2164,9 +2161,9 @@ version = "0.21.71" [[deps.VortexStepMethod]] deps = ["Colors", "DefaultApplication", "DelimitedFiles", "DifferentiationInterface", "FiniteDiff", "Interpolations", "LaTeXStrings", "LinearAlgebra", "Logging", "Measures", "NonlinearSolve", "Parameters", "Pkg", "PreallocationTools", "PrecompileTools", "SciMLBase", "Serialization", "StaticArrays", "Statistics", "Timers", "Xfoil", "YAML"] -git-tree-sha1 = "fd7616ebd3dc9fc2c0e5e6abf7ba641aa52be444" +git-tree-sha1 = "a0fa008357206434acfa2acc588b84764478df62" uuid = "ed3cd733-9f0f-46a9-93e0-89b8d4998dd9" -version = "1.2.5" +version = "1.2.6" [deps.VortexStepMethod.extensions] VortexStepMethodControlPlotsExt = "ControlPlots" @@ -2216,9 +2213,9 @@ version = "0.5.0" [[deps.YAML]] deps = ["Base64", "Dates", "Printf", "StringEncodings"] -git-tree-sha1 = "b46894beba6c05cd185d174654479aaec09ea6b1" +git-tree-sha1 = "2f58ac39f64b41fb812340347525be3b590cce3b" uuid = "ddb6d928-2868-570f-bddf-ab3f9cf99eb6" -version = "0.4.13" +version = "0.4.14" [[deps.Zlib_jll]] deps = ["Libdl"] diff --git a/Manifest-v1.11.toml.default.license b/Manifest-v1.11.toml.default.license deleted file mode 100644 index c52bffc3b..000000000 --- a/Manifest-v1.11.toml.default.license +++ /dev/null @@ -1,2 +0,0 @@ -SPDX-FileCopyrightText: 2025 Uwe Fechner -SPDX-License-Identifier: MIT diff --git a/Manifest-v1.11.toml.license b/Manifest-v1.11.toml.license deleted file mode 100644 index fa17de081..000000000 --- a/Manifest-v1.11.toml.license +++ /dev/null @@ -1,3 +0,0 @@ -SPDX-FileCopyrightText: 2025 Uwe Fechner - -SPDX-License-Identifier: MIT diff --git a/Project.toml b/Project.toml index 6354c8575..0cf61b1cc 100644 --- a/Project.toml +++ b/Project.toml @@ -61,15 +61,15 @@ JLD2 = "0.5" KitePodModels = "0.3.8" KiteUtils = "0.10.5" LaTeXStrings = "1.4.0" -LinearSolve = "~2.39.0" LinearAlgebra = "1.10, 1.11" -ModelingToolkit = "~9.72.0" +LinearSolve = "~2.39.0" +ModelingToolkit = "~9.78.0" NLSolversBase = "~7.8.3" NLsolve = "4.5" -NonlinearSolve = "=4.5.0" -OrdinaryDiffEqBDF = "1.4.0" -OrdinaryDiffEqCore = "1.22.0" -OrdinaryDiffEqNonlinearSolve = "1.6.0" +NonlinearSolve = "4.8.0" +OrdinaryDiffEqBDF = "1.5.0" +OrdinaryDiffEqCore = "1.23.0" +OrdinaryDiffEqNonlinearSolve = "1.6.1" PackageCompiler = "2.0" Parameters = "0.12" Pkg = "1.10, 1.11" @@ -85,10 +85,10 @@ Statistics = "1" StatsBase = "0.34" Sundials = "4.24" SymbolicIndexingInterface = "0.3" -SymbolicUtils = "~3.25" +SymbolicUtils = "3.25" Test = "1" Timers = "0.1.5" -VortexStepMethod = "1.2.5" +VortexStepMethod = "1.2.6" WinchModels = "0.3.6" julia = "1.10, 1.11" @@ -97,8 +97,8 @@ Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" Colors = "5ae59095-9a9b-59fe-a467-6f913c188581" ControlPlots = "23c2ee80-7a9e-4350-b264-8e670f12517c" -Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2" +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" PackageCompiler = "9b87118b-4619-50d2-8e1e-99f35a4d4d9d" diff --git a/REUSE.toml b/REUSE.toml new file mode 100644 index 000000000..9ac0b1e5b --- /dev/null +++ b/REUSE.toml @@ -0,0 +1,43 @@ +version = 1 + +[[annotations]] +path = "data/prob*" + +SPDX-FileCopyrightText = "2025 Bart van de Lint" +SPDX-License-Identifier = "MPL-2.0" + +[[annotations]] +path = "data/ram_air_kite*" + +SPDX-FileCopyrightText = "2025 Bart van de Lint" +SPDX-License-Identifier = "MPL-2.0" + +[[annotations]] +path = "Manifest*" + +SPDX-FileCopyrightText = "2025 Uwe Fechner, Bart van de Lint" +SPDX-License-Identifier = "MIT" + +[[annotations]] +path = "Project.toml" + +SPDX-FileCopyrightText = "2025 Uwe Fechner, Bart van de Lint" +SPDX-License-Identifier = "MIT" + +[[annotations]] +path = "CITATION.cff" + +SPDX-FileCopyrightText = "2025 Uwe Fechner" +SPDX-License-Identifier = "MIT" + +[[annotations]] +path = ".gitignore" + +SPDX-FileCopyrightText = "2025 Uwe Fechner, Bart van de Lint" +SPDX-License-Identifier = "MIT" + +[[annotations]] +path = "docs/src/*.png" + +SPDX-FileCopyrightText = "2025 Uwe Fechner" +SPDX-License-Identifier = "CC-BY-4.0" diff --git a/bin/create_sys_image b/bin/create_sys_image index b6ca7b53b..af07c20cc 100755 --- a/bin/create_sys_image +++ b/bin/create_sys_image @@ -1,16 +1,42 @@ #!/bin/bash -eu +# SPDX-FileCopyrightText: 2025 Uwe Fechner +# SPDX-License-Identifier: MIT update=false -if [[ $# -gt 0 ]]; then - if [[ $1 != "--update" ]]; then - echo "Invalid parameter! Use:" - echo "./create_sys_image" - echo "or" - echo "./create_sys_image --update" - exit 1 - else - update=true - fi +julia_version="" + +# Parse command line arguments +while [[ $# -gt 0 ]]; do + case $1 in + --update) + update=true + shift + ;; + --version) + if [[ $# -lt 2 ]]; then + echo "Error: --version requires a version argument" + exit 1 + fi + julia_version="$2" + shift 2 + ;; + *) + echo "Invalid parameter! Use:" + echo "./create_sys_image" + echo "or" + echo "./create_sys_image --update" + echo "or" + echo "./create_sys_image --version 1.x" + exit 1 + ;; + esac +done + +# Remove the version switching code +if [[ -n "$julia_version" ]]; then + JULIA_CMD="julia +$julia_version" +else + JULIA_CMD="julia" fi if [[ $(basename $(pwd)) == "bin" ]]; then @@ -68,7 +94,7 @@ if [[ $totalmem -lt 27000 ]]; then export JULIA_IMAGE_THREADS=1 fi -julia_version=$(julia --version | awk '{print($3)}') +julia_version=$($JULIA_CMD --version | awk '{print($3)}') julia_major=${julia_version:0:3} if [[ $julia_major == "1.1" ]]; then julia_major=${julia_version:0:4} @@ -81,7 +107,7 @@ fi if test -f "kps-image-${julia_major}-${branch}.so"; then mv bin/kps-image-${julia_major}-${branch}.so kps-image-${julia_major}-${branch}.so.bak fi -julia --startup-file=no -e "using Pkg; Pkg.add(\"TestEnv\")" +$JULIA_CMD --startup-file=no -e "using Pkg; Pkg.add(\"TestEnv\")" if [[ $update == true ]]; then echo "Updating packages..." @@ -98,11 +124,11 @@ if [[ $update == true ]]; then mv Manifest-v1.11.toml Manifest-v1.11.toml.bak fi fi - julia --startup-file=no --pkgimages=no --project -e "using TestEnv; TestEnv.activate(); using Pkg; Pkg.add(\"PyCall\"); Pkg.build(\"PyCall\")" + $JULIA_CMD --startup-file=no --pkgimages=no --project -e "using TestEnv; TestEnv.activate(); using Pkg; Pkg.add(\"PyCall\"); Pkg.build(\"PyCall\")" if [[ $PYTHON == "" ]]; then - julia --startup-file=no --pkgimages=no --project -e "using TestEnv; TestEnv.activate(); using Pkg; Pkg.add(\"Conda\"); using Conda; Conda.add(\"matplotlib\"); using ControlPlots" + $JULIA_CMD --startup-file=no --pkgimages=no --project -e "using TestEnv; TestEnv.activate(); using Pkg; Pkg.add(\"Conda\"); using Conda; Conda.add(\"matplotlib\"); using ControlPlots" fi - julia --startup-file=no --pkgimages=no --project -e "using Pkg; Pkg.update()" + $JULIA_CMD --startup-file=no --pkgimages=no --project -e "using Pkg; Pkg.update()" if [[ $julia_major == "1.10" ]]; then mv Manifest.toml Manifest-v1.10.toml else @@ -118,10 +144,10 @@ else fi fi if [[ $update == true ]]; then - julia --startup-file=no --pkgimages=no --project -e "include(\"./test/create_sys_image.jl\");" + $JULIA_CMD --startup-file=no --pkgimages=no --project -e "include(\"./test/create_sys_image.jl\");" else - julia --startup-file=no --pkgimages=no --project -e "using Pkg; Pkg.instantiate();" - julia --startup-file=no --pkgimages=no --project -e "using TestEnv; TestEnv.activate(); include(\"./test/create_sys_image.jl\");" + $JULIA_CMD --startup-file=no --pkgimages=no --project -e "using Pkg; Pkg.instantiate();" + $JULIA_CMD --startup-file=no --pkgimages=no --project -e "using TestEnv; TestEnv.activate(); include(\"./test/create_sys_image.jl\");" fi if [[ $branch != "" ]]; then @@ -133,19 +159,20 @@ if test -f $SOFILE; then mv $SOFILE $SOFILE.bak fi mv kps-image_tmp.so $SOFILE -julia --startup-file=no --project -e "using Pkg; Pkg.precompile(); Pkg.resolve()" +$JULIA_CMD --startup-file=no --project -e "using Pkg; Pkg.precompile(); Pkg.resolve()" if [ -d src ]; then cd src touch *.jl # make sure all modules get recompiled in the next step cd .. fi -rm data/*.bin -rm data/*.bin.default +rm -f data/*.bin +rm -f data/*.bin.default echo "Precompiling package KiteModels..." -julia --startup-file=no --project -J $SOFILE -e "using KiteModels, KitePodModels, KiteUtils" +$JULIA_CMD --startup-file=no --project -J $SOFILE -e "using KiteModels, KitePodModels, KiteUtils" echo "Precompiling package ControlPlots..." if [[ $branch != "" ]]; then - julia --startup-file=no --project -J $SOFILE -e "using TestEnv; TestEnv.activate(); using ControlPlots; using VortexStepMethod" + $JULIA_CMD --startup-file=no --project -J $SOFILE -e "using TestEnv; TestEnv.activate(); using ControlPlots; using VortexStepMethod" else - julia --startup-file=no --project -J $SOFILE -e "using ControlPlots; using VortexStepMethod" + $JULIA_CMD --startup-file=no --project -J $SOFILE -e "using ControlPlots; using VortexStepMethod" fi + diff --git a/bin/create_sys_image.license b/bin/create_sys_image.license deleted file mode 100644 index fa17de081..000000000 --- a/bin/create_sys_image.license +++ /dev/null @@ -1,3 +0,0 @@ -SPDX-FileCopyrightText: 2025 Uwe Fechner - -SPDX-License-Identifier: MIT diff --git a/bin/create_xz_file b/bin/create_xz_file index d788b6edc..d8930dbe9 100755 --- a/bin/create_xz_file +++ b/bin/create_xz_file @@ -1,4 +1,6 @@ #!/bin/bash -eu +# SPDX-FileCopyrightText: 2025 Uwe Fechner +# SPDX-License-Identifier: MIT if [[ $(basename $(pwd)) == "bin" ]]; then cd .. @@ -22,18 +24,24 @@ fi echo "Creating xz file for prob_dynamic_v${julia_major}_3_seg.bin.default.xz ..." if [[ $julia_major == "1.10" ]]; then - rm -f data/prob_dynamic_1.10_3_seg.bin.* - echo "using KiteModels..." + rm -f data/prob_1.10_ram_dynamic_3_seg.bin.default + rm -f data/prob_1.10_ram_dynamic_3_seg.bin.default.xz + echo "using KiteModels. Please ignore the error message about the missing input file!" printf ' ' >> src/precompile.jl - julia --project -J $SOFILE -e 'using KiteModels' - cp data/prob_dynamic_1.10_3_seg.bin data/prob_dynamic_1.10_3_seg.bin.default - xz data/prob_dynamic_1.10_3_seg.bin.default + julia --project -J $SOFILE -e 'using KiteModels; include("test/create_xz_file.jl")' + cp data/prob_1.10_ram_dynamic_3_seg.bin data/prob_1.10_ram_dynamic_3_seg.bin.default + xz data/prob_1.10_ram_dynamic_3_seg.bin.default else - rm -f data/prob_dynamic_1.11_3_seg.bin.default.xz - cp data/prob_dynamic_1.11_3_seg.bin data/prob_dynamic_1.11_3_seg.bin.default - xz data/prob_dynamic_1.11_3_seg.bin.default + rm -f data/prob_1.11_ram_dynamic_3_seg.bin.default + rm -f data/prob_1.11_ram_dynamic_3_seg.bin.default.xz + echo "using KiteModels. Please ignore the error message about the missing input file!" + printf ' ' >> src/precompile.jl + julia --project -J $SOFILE -e 'using KiteModels; include("test/create_xz_file.jl")' + cp data/prob_1.11_ram_dynamic_3_seg.bin data/prob_1.11_ram_dynamic_3_seg.bin.default + xz data/prob_1.11_ram_dynamic_3_seg.bin.default fi echo "using KiteModels..." truncate -s -1 src/precompile.jl julia --project -J $SOFILE -e 'using KiteModels' +echo "Done creating prob_dynamic_v${julia_major}_3_seg.bin.default.xz !" diff --git a/bin/pull b/bin/pull index 6c49cafa0..d63af7b8b 100755 --- a/bin/pull +++ b/bin/pull @@ -1,4 +1,6 @@ -#!/bin/bash +#!/bin/bash -eu +# SPDX-FileCopyrightText: 2025 Uwe Fechner +# SPDX-License-Identifier: MIT if [[ $(basename $(pwd)) == "bin" ]]; then cd .. diff --git a/bin/pull.license b/bin/pull.license deleted file mode 100644 index fa17de081..000000000 --- a/bin/pull.license +++ /dev/null @@ -1,3 +0,0 @@ -SPDX-FileCopyrightText: 2025 Uwe Fechner - -SPDX-License-Identifier: MIT diff --git a/bin/run_julia b/bin/run_julia index 205f63b7e..4a70dad78 100755 --- a/bin/run_julia +++ b/bin/run_julia @@ -1,4 +1,6 @@ -#!/bin/bash +#!/bin/bash -eu +# SPDX-FileCopyrightText: 2025 Uwe Fechner +# SPDX-License-Identifier: MIT if [[ $(basename $(pwd)) == "bin" ]]; then cd .. @@ -29,15 +31,15 @@ fi if [[ $branch != "" ]]; then if test -f "bin/kps-image-${julia_major}-${branch}.so"; then echo "Found system image!" - julia -J bin/kps-image-${julia_major}-${branch}.so -t 1 $GCT --project -i -e 'using KiteModels' + julia +${julia_major} -J bin/kps-image-${julia_major}-${branch}.so -t 1 $GCT --project -i -e 'using KiteModels' else - julia $GCT --project -i # -e 'using KiteModels' + julia +${julia_major} $GCT --project -i # -e 'using KiteModels' fi else if test -f "bin/kps-image-${julia_major}.so"; then echo "Found system image!" - julia -J bin/kps-image-${julia_major}.so -t 1 $GCT --project -i -e 'using KiteModels' + julia +${julia_major} -J bin/kps-image-${julia_major}.so -t 1 $GCT --project -i -e 'using KiteModels' else - julia $GCT --project -i # -e 'using KiteModels' + julia +${julia_major} $GCT --project -i # -e 'using KiteModels' fi -fi +fi \ No newline at end of file diff --git a/bin/run_julia.license b/bin/run_julia.license deleted file mode 100644 index fa17de081..000000000 --- a/bin/run_julia.license +++ /dev/null @@ -1,3 +0,0 @@ -SPDX-FileCopyrightText: 2025 Uwe Fechner - -SPDX-License-Identifier: MIT diff --git a/bin/update_default_manifest b/bin/update_default_manifest index 95318cdde..dcac8f3d2 100755 --- a/bin/update_default_manifest +++ b/bin/update_default_manifest @@ -1,4 +1,6 @@ #!/bin/bash -eu +# SPDX-FileCopyrightText: 2025 Uwe Fechner +# SPDX-License-Identifier: MIT if [[ $(basename $(pwd)) == "bin" ]]; then cd .. diff --git a/bin/update_default_manifest.license b/bin/update_default_manifest.license deleted file mode 100644 index fa17de081..000000000 --- a/bin/update_default_manifest.license +++ /dev/null @@ -1,3 +0,0 @@ -SPDX-FileCopyrightText: 2025 Uwe Fechner - -SPDX-License-Identifier: MIT diff --git a/data/prob_1.10_ram_dynamic_3_seg.bin.default.xz b/data/prob_1.10_ram_dynamic_3_seg.bin.default.xz new file mode 100644 index 000000000..f64ae4cbd Binary files /dev/null and b/data/prob_1.10_ram_dynamic_3_seg.bin.default.xz differ diff --git a/data/prob_1.11_ram_dynamic_3_seg.bin.default.xz b/data/prob_1.11_ram_dynamic_3_seg.bin.default.xz new file mode 100644 index 000000000..edae9ec8f Binary files /dev/null and b/data/prob_1.11_ram_dynamic_3_seg.bin.default.xz differ diff --git a/data/prob_dynamic_1.10_3_seg.bin.default.xz.license b/data/prob_dynamic_1.10_3_seg.bin.default.xz.license deleted file mode 100644 index 74717d8f7..000000000 --- a/data/prob_dynamic_1.10_3_seg.bin.default.xz.license +++ /dev/null @@ -1,3 +0,0 @@ -SPDX-FileCopyrightText: 2025 Bart van de Lint - -SPDX-License-Identifier: CC-BY-4.0 diff --git a/data/prob_dynamic_1.11_3_seg.bin.default.license b/data/prob_dynamic_1.11_3_seg.bin.default.license deleted file mode 100644 index 74717d8f7..000000000 --- a/data/prob_dynamic_1.11_3_seg.bin.default.license +++ /dev/null @@ -1,3 +0,0 @@ -SPDX-FileCopyrightText: 2025 Bart van de Lint - -SPDX-License-Identifier: CC-BY-4.0 diff --git a/data/prob_dynamic_1.11_3_seg.bin.default.xz b/data/prob_dynamic_1.11_3_seg.bin.default.xz deleted file mode 100644 index 477f897e4..000000000 Binary files a/data/prob_dynamic_1.11_3_seg.bin.default.xz and /dev/null differ diff --git a/data/prob_dynamic_1.11_3_seg.bin.default.xz.license b/data/prob_dynamic_1.11_3_seg.bin.default.xz.license deleted file mode 100644 index 74717d8f7..000000000 --- a/data/prob_dynamic_1.11_3_seg.bin.default.xz.license +++ /dev/null @@ -1,3 +0,0 @@ -SPDX-FileCopyrightText: 2025 Bart van de Lint - -SPDX-License-Identifier: CC-BY-4.0 diff --git a/data/ram_air_kite_body.obj.license b/data/ram_air_kite_body.obj.license deleted file mode 100644 index 74717d8f7..000000000 --- a/data/ram_air_kite_body.obj.license +++ /dev/null @@ -1,3 +0,0 @@ -SPDX-FileCopyrightText: 2025 Bart van de Lint - -SPDX-License-Identifier: CC-BY-4.0 diff --git a/data/ram_air_kite_foil.dat.license b/data/ram_air_kite_foil.dat.license deleted file mode 100644 index 74717d8f7..000000000 --- a/data/ram_air_kite_foil.dat.license +++ /dev/null @@ -1,3 +0,0 @@ -SPDX-FileCopyrightText: 2025 Bart van de Lint - -SPDX-License-Identifier: CC-BY-4.0 diff --git a/data/ram_air_kite_foil_cd_polar.csv.license b/data/ram_air_kite_foil_cd_polar.csv.license deleted file mode 100644 index 74717d8f7..000000000 --- a/data/ram_air_kite_foil_cd_polar.csv.license +++ /dev/null @@ -1,3 +0,0 @@ -SPDX-FileCopyrightText: 2025 Bart van de Lint - -SPDX-License-Identifier: CC-BY-4.0 diff --git a/data/ram_air_kite_foil_cl_polar.csv.license b/data/ram_air_kite_foil_cl_polar.csv.license deleted file mode 100644 index 74717d8f7..000000000 --- a/data/ram_air_kite_foil_cl_polar.csv.license +++ /dev/null @@ -1,3 +0,0 @@ -SPDX-FileCopyrightText: 2025 Bart van de Lint - -SPDX-License-Identifier: CC-BY-4.0 diff --git a/data/ram_air_kite_foil_cm_polar.csv.license b/data/ram_air_kite_foil_cm_polar.csv.license deleted file mode 100644 index 74717d8f7..000000000 --- a/data/ram_air_kite_foil_cm_polar.csv.license +++ /dev/null @@ -1,3 +0,0 @@ -SPDX-FileCopyrightText: 2025 Bart van de Lint - -SPDX-License-Identifier: CC-BY-4.0 diff --git a/data/settings_ram.yaml b/data/settings_ram.yaml index 6883631eb..9c4051765 100644 --- a/data/settings_ram.yaml +++ b/data/settings_ram.yaml @@ -15,7 +15,7 @@ initial: v_reel_out: 0.0 # initial reel out speed [m/s] solver: - abs_tol: 0.006 # absolute tolerance of the DAE solver [m, m/s] + abs_tol: 0.01 # absolute tolerance of the DAE solver [m, m/s] rel_tol: 0.01 # relative tolerance of the DAE solver [-] solver: "DFBDF" # DAE solver, IDA or DImplicitEuler, DFBDF linear_solver: "GMRES" # can be GMRES or Dense or LapackDense (only for IDA) @@ -25,7 +25,7 @@ solver: kite: model: "data/ram_air_kite_body.obj" # 3D model of the kite foil_file: "data/ram_air_kite_foil.dat" # filename for the foil shape - physical_model: "RamAirKite" # name of the kite model to use (KPS3, KPS4 or RamAirKite) + physical_model: "ram" # name of the kite model to use (KPS3, KPS4 or RamAirKite) top_bridle_points: # top bridle points that are not on the kite body in CAD frame - [0.290199, 0.784697, -2.61305] - [0.392683, 0.785271, -2.61201] diff --git a/docs/src/4-point-kite.png.license b/docs/src/4-point-kite.png.license deleted file mode 100644 index 48d71b826..000000000 --- a/docs/src/4-point-kite.png.license +++ /dev/null @@ -1,3 +0,0 @@ -SPDX-FileCopyrightText: 2025 Uwe Fechner - -SPDX-License-Identifier: CC-BY-4.0 diff --git a/docs/src/initial_state.png.license b/docs/src/initial_state.png.license deleted file mode 100644 index 48d71b826..000000000 --- a/docs/src/initial_state.png.license +++ /dev/null @@ -1,3 +0,0 @@ -SPDX-FileCopyrightText: 2025 Uwe Fechner - -SPDX-License-Identifier: CC-BY-4.0 diff --git a/docs/src/initial_state_4p.png.license b/docs/src/initial_state_4p.png.license deleted file mode 100644 index 48d71b826..000000000 --- a/docs/src/initial_state_4p.png.license +++ /dev/null @@ -1,3 +0,0 @@ -SPDX-FileCopyrightText: 2025 Uwe Fechner - -SPDX-License-Identifier: CC-BY-4.0 diff --git a/docs/src/kite.png.license b/docs/src/kite.png.license deleted file mode 100644 index 48d71b826..000000000 --- a/docs/src/kite.png.license +++ /dev/null @@ -1,3 +0,0 @@ -SPDX-FileCopyrightText: 2025 Uwe Fechner - -SPDX-License-Identifier: CC-BY-4.0 diff --git a/docs/src/kite_power_tools.png.license b/docs/src/kite_power_tools.png.license deleted file mode 100644 index 48d71b826..000000000 --- a/docs/src/kite_power_tools.png.license +++ /dev/null @@ -1,3 +0,0 @@ -SPDX-FileCopyrightText: 2025 Uwe Fechner - -SPDX-License-Identifier: CC-BY-4.0 diff --git a/docs/src/kps4.png.license b/docs/src/kps4.png.license deleted file mode 100644 index 48d71b826..000000000 --- a/docs/src/kps4.png.license +++ /dev/null @@ -1,3 +0,0 @@ -SPDX-FileCopyrightText: 2025 Uwe Fechner - -SPDX-License-Identifier: CC-BY-4.0 diff --git a/docs/src/kps4_hires.png.license b/docs/src/kps4_hires.png.license deleted file mode 100644 index 48d71b826..000000000 --- a/docs/src/kps4_hires.png.license +++ /dev/null @@ -1,3 +0,0 @@ -SPDX-FileCopyrightText: 2025 Uwe Fechner - -SPDX-License-Identifier: CC-BY-4.0 diff --git a/docs/src/reelout_force_1p.png.license b/docs/src/reelout_force_1p.png.license deleted file mode 100644 index 48d71b826..000000000 --- a/docs/src/reelout_force_1p.png.license +++ /dev/null @@ -1,3 +0,0 @@ -SPDX-FileCopyrightText: 2025 Uwe Fechner - -SPDX-License-Identifier: CC-BY-4.0 diff --git a/docs/src/reelout_force_4p.png.license b/docs/src/reelout_force_4p.png.license deleted file mode 100644 index 48d71b826..000000000 --- a/docs/src/reelout_force_4p.png.license +++ /dev/null @@ -1,3 +0,0 @@ -SPDX-FileCopyrightText: 2025 Uwe Fechner - -SPDX-License-Identifier: CC-BY-4.0 diff --git a/docs/src/test_plan.md b/docs/src/test_plan.md index ea7e7c436..c4caa748f 100644 --- a/docs/src/test_plan.md +++ b/docs/src/test_plan.md @@ -2,20 +2,27 @@ On AMD Ryzen 7 7840U. ### Pass criteria for the performance: +First run: +- simplifying the system < 45s + +Second run: - system initialized < 15s - total time without plotting < 35s ### Develop branch +In Bash: ``` juliaup default 1.11 git checkout develop cd bin -./create_sys_image --update +./update_default_manifest +./update_xz_file cd .. -rm data/prob_1.11*.bin +git commit -m "Update default manifest" Manifest-v1.11.toml.default +git push jl ``` -On a laptop using +Now in Julia: ``` include("examples/ram_air_kite.jl") ``` @@ -23,72 +30,23 @@ Expected output: ``` julia> include("examples/ram_air_kite.jl") [ Info: Loading packages -Time elapsed: 2.843903577 s +Time elapsed: 2.863292943 s [ Info: Creating wing, aero, vsm_solver, point_system and s: -Time elapsed: 11.780494151 s -[ Info: Creating ODESystem - 4.189182 seconds (7.83 M allocations: 200.190 MiB, 0.69% gc time, 24.40% compilation time: 6% of which was recompilation) -[ Info: Simplifying the system - 38.041611 seconds (304.71 M allocations: 10.010 GiB, 2.12% gc time, 25.50% compilation time: 18% of which was recompilation) -[ Info: Creating ODEProblem - 84.739970 seconds (640.11 M allocations: 20.794 GiB, 1.83% gc time, 22.66% compilation time: 12% of which was recompilation) -[ Info: Initialized integrator in 22.253394034 seconds +Time elapsed: 4.51880941 s +[ Info: Initialized integrator in 6.000656365 seconds [ Info: System initialized at: -Time elapsed: 193.93390819 s +Time elapsed: 13.392478822 s [ Info: Total time without plotting: -Time elapsed: 211.095669066 s +Time elapsed: 30.847251362 s ┌ Info: Performance: -│ times_realtime = 8.017735204615859 -└ integrator_times_realtime = 60.411279264296 +│ times_realtime = 8.426564704784319 +└ integrator_times_realtime = 62.47122497288717 ``` -### Main branch +### NLSolve test ``` -julia> include("examples/ram_air_kite.jl") -[ Info: Loading packages -Time elapsed: 2.9059458 s -[ Info: Creating wing, aero, vsm_solver, point_system and s: -Time elapsed: 11.790243385 s -┌ Warning: Failure to deserialize /home/ufechner/repos/KiteModels.jl/data/prob_dynamic_1.11_3_seg.bin ! -└ @ KiteModels ~/repos/KiteModels.jl/src/ram_air_kite.jl:519 -[ Info: Rebuilding the system. This can take some minutes... -[ Info: Creating ODESystem - 4.775993 seconds (8.31 M allocations: 201.174 MiB, 0.64% gc time, 20.63% compilation time: 6% of which was recompilation) -[ Info: Simplifying the system - 43.592286 seconds (329.71 M allocations: 10.926 GiB, 2.06% gc time, 25.35% compilation time: 18% of which was recompilation) -[ Info: Creating ODEProblem - 95.744396 seconds (682.85 M allocations: 23.818 GiB, 2.08% gc time, 17.67% compilation time: 14% of which was recompilation) -[ Info: Initialized integrator in 23.951264242 seconds -[ Info: System initialized at: -Time elapsed: 212.672152747 s -[ Info: Total time without plotting: -Time elapsed: 232.546116606 s -┌ Info: Performance: -│ times_realtime = 4.415094072484458 -└ integrator_times_realtime = 14.752422253318498 +include("mwes/mwe_26.jl") +include("mwes/mwe_26.jl") +include("mwes/mwe_26.jl") ``` - -### Update precompiled code -``` -cd bin -./update_default_manifest -cd .. -cd data -cp prob_dynamic_1.11_3_seg.bin prob_dynamic_1.11_3_seg.bin.default -cd .. -echo " " >> src/precompile.jl -``` -## Test again -julia> include("examples/ram_air_kite.jl") -[ Info: Loading packages -Time elapsed: 2.807525651 s -[ Info: Creating wing, aero, vsm_solver, point_system and s: -Time elapsed: 4.499499709 s -[ Info: Initialized integrator in 5.94955637 seconds -[ Info: System initialized at: -Time elapsed: 13.388212381 s -[ Info: Total time without plotting: -Time elapsed: 33.223129728 s -┌ Info: Performance: -│ times_realtime = 4.4603427235287265 -└ integrator_times_realtime = 15.033714698785786 \ No newline at end of file +Each time the force, that is printed must be > 9000 N. diff --git a/docs/src/winch.png.license b/docs/src/winch.png.license deleted file mode 100644 index 48d71b826..000000000 --- a/docs/src/winch.png.license +++ /dev/null @@ -1,3 +0,0 @@ -SPDX-FileCopyrightText: 2025 Uwe Fechner - -SPDX-License-Identifier: CC-BY-4.0 diff --git a/examples/discrete_jacobian.jl b/examples/discrete_jacobian.jl new file mode 100644 index 000000000..68eb98e50 --- /dev/null +++ b/examples/discrete_jacobian.jl @@ -0,0 +1,162 @@ +# Copyright (c) 2025 Bart van de Lint +# SPDX-License-Identifier: MPL-2.0 + +#= +This example analyzes the input-output behavior of a ram air kite model by: + +1. Creating a ram air kite model and initializing it at 60° elevation +2. Stabilizing the system +3. Testing the system response to different control inputs: + - Left steering line torque + - Right steering line torque + +The script plots relationships between steering inputs and resulting angular velocities, +providing insight into the kite's steering behavior and control characteristics. +=# + +using Timers +tic() +@info "Loading packages " + +using KiteModels, LinearAlgebra, Statistics, OrdinaryDiffEqCore, OrdinaryDiffEqNonlinearSolve, OrdinaryDiffEqBDF +using DifferentiationInterface, FiniteDiff, SparseArrays, ADTypes +using ModelingToolkit: setu, getu, setp, getp +using ModelingToolkit + +NOISE = 0.1 +PLOT = true +if PLOT + using Pkg + if ! ("LaTeXStrings" ∈ keys(Pkg.project().dependencies)) + using TestEnv; TestEnv.activate() + end + using ControlPlots, LaTeXStrings + import ControlPlots: plot +end +toc() + + +include(joinpath(@__DIR__, "plotting.jl")) + +# Simulation parameters +dt = 0.05 +total_time = 10 # Longer simulation to see oscillations +vsm_interval = 3 +steps = Int(round(total_time / dt)) + +# Steering parameters +steering_freq = 1/2 # Hz - full left-right cycle frequency +steering_magnitude = 10.0 # Magnitude of steering input [Nm] + +# Initialize model +set = load_settings("system_ram.yaml") +set.segments = 2 +set_values = [-50, 0.0, 0.0] # Set values of the torques of the three winches. [Nm] +set.quasi_static = true +set.physical_model = "ram" + +@info "Creating wing, aero, vsm_solver, point_system and s:" +s = RamAirKite(set) +s.set.abs_tol = 1e-3 +s.set.rel_tol = 1e-3 +toc() + +# init_Q_b_w, R_b_w = KiteModels.measure_to_q(measure) +# init_kite_pos = init!(s.point_system, s.set, R_b_w) +# plot(s.point_system, 0.0; zoom=false, front=true) + +measure = Measurement() + +# Initialize at elevation +measure.sphere_pos .= deg2rad.([70.0 70.0; 1.0 -1.0]) +KiteModels.init_sim!(s, measure; + adaptive=false, remake=false, reload=true, + solver=FBDF(nlsolve=OrdinaryDiffEqNonlinearSolve.NLNewton(relax=0.8, max_iter=1000)) +) +OrdinaryDiffEqCore.set_proposed_dt!(s.integrator, dt) +sys = s.sys + +@info "System initialized at:" +toc() + +sys = s.sys +# # Stabilize system +s.integrator.ps[sys.stabilize] = true +for i in 1:10÷dt + next_step!(s; dt, vsm_interval=1) +end +s.integrator.ps[sys.stabilize] = false +# plot(s, 0.0; zoom=true, front=false) + +# Function to step simulation with input u +function f(x, u, _, p) + (s, set_x, set_u, get_x, dt) = p + set_x(s.integrator, x) + set_u(s.integrator, u) + OrdinaryDiffEqCore.reinit!(s.integrator, s.integrator.u; reinit_dae=false) + OrdinaryDiffEqCore.step!(s.integrator, dt) + xnext = get_x(s.integrator) + @show norm(xnext) + return xnext +end + +# Get initial state +x_vec = KiteModels.get_unknowns(s) +set_ix = setu(s.integrator, Initial.(x_vec)) # set_ix might not be needed anymore if prob is removed +set_x = setu(s.integrator, x_vec) +set_u = setu(s.integrator, collect(sys.set_values)) +get_u = getu(s.integrator, collect(sys.set_values)) +get_x = getu(s.integrator, x_vec) +x0 = get_x(s.integrator) +u0 = get_u(s.integrator) +p = (s, set_x, set_u, get_x, dt) + +@show x0 +x0 .+= (rand(length(x0)) .- 0.5) * NOISE +set_ix(s.integrator, x0) +@time OrdinaryDiffEqCore.reinit!(s.integrator) +@time OrdinaryDiffEqCore.reinit!(s.integrator) +@time OrdinaryDiffEqCore.reinit!(s.integrator) + +# use_scc=false: +# 0.026057 seconds (2.04 k allocations: 712.352 KiB) +# use_ssc=true: +# 0.017217 seconds (4.08 k allocations: 574.000 KiB) + +f_x(x) = f(x, u0, nothing, p) +f_u(u) = f(x0, u, nothing, p) + +x_idxs = Dict{Num, Int}() +for (idx, sym) in enumerate(x_vec) + x_idxs[sym] = idx +end + +m, n = length(x_vec), length(x_vec) +active_x = [ + [x_idxs[sys.free_twist_angle[i]] for i in 1:4] + [x_idxs[sys.ω_b[i]] for i in 1:3] + [x_idxs[sys.kite_pos[i]] for i in 1:3] + [x_idxs[sys.kite_vel[i]] for i in 1:3] +] +# Collect all (row, col) pairs +rows = Int[] +cols = Int[] +for i in active_x, j in active_x + push!(rows, i) + push!(cols, j) +end +# Build the sparse boolean matrix +S = sparse(rows, cols, trues(length(rows)), m, n) + + +backend = AutoFiniteDiff() +sparse_backend = AutoSparse( + backend; + sparsity_detector=ADTypes.KnownJacobianSparsityDetector(S) +) +# jac_prep = prepare_jacobian(f_x, backend, x0) +# @time J = jacobian(f_x, jac_prep, backend, x0) +# @time J = jacobian(f_x, jac_prep, backend, x0) +# @show norm(J) +# display(J) + diff --git a/examples/input_output_function.jl b/examples/input_output_function.jl new file mode 100644 index 000000000..7a46875d9 --- /dev/null +++ b/examples/input_output_function.jl @@ -0,0 +1,211 @@ +# Copyright (c) 2025 Bart van de Lint +# SPDX-License-Identifier: MPL-2.0 + +#= +This example analyzes the input-output behavior of a ram air kite model by: + +1. Creating a ram air kite model and initializing it at 60° elevation +2. Stabilizing the system +3. Testing the system response to different control inputs: + - Left steering line torque + - Right steering line torque + +The script plots relationships between steering inputs and resulting angular velocities, +providing insight into the kite's steering behavior and control characteristics. +=# + +using Timers +tic() +@info "Loading packages " + +using KiteModels, LinearAlgebra, Statistics, OrdinaryDiffEqCore, OrdinaryDiffEqNonlinearSolve, OrdinaryDiffEqBDF +using ModelingToolkit: setu, getu, setp, getp +using ModelingToolkit + +PLOT = true +if PLOT + using Pkg + if ! ("LaTeXStrings" ∈ keys(Pkg.project().dependencies)) + using TestEnv; TestEnv.activate() + end + using ControlPlots, LaTeXStrings + import ControlPlots: plot +end +toc() + + +include(joinpath(@__DIR__, "plotting.jl")) + +# Simulation parameters +dt = 0.05 +total_time = 10 # Longer simulation to see oscillations +vsm_interval = 3 +steps = Int(round(total_time / dt)) + +# Steering parameters +steering_freq = 1/2 # Hz - full left-right cycle frequency +steering_magnitude = 10.0 # Magnitude of steering input [Nm] + +# Initialize model +set = load_settings("system_ram.yaml") +set.segments = 2 +set_values = [-50, 0.0, 0.0] # Set values of the torques of the three winches. [Nm] +set.quasi_static = true +set.physical_model = "ram" + +@info "Creating wing, aero, vsm_solver, point_system and s:" +s = RamAirKite(set) +s.set.abs_tol = 1e-2 +s.set.rel_tol = 1e-2 +toc() + +# init_Q_b_w, R_b_w = KiteModels.measure_to_q(measure) +# init_kite_pos = init!(s.point_system, s.set, R_b_w) +# plot(s.point_system, 0.0; zoom=false, front=true) + +measure = Measurement() + +# Initialize at elevation +measure.sphere_pos .= deg2rad.([70.0 70.0; 1.0 -1.0]) +KiteModels.init_sim!(s, measure; + adaptive=false, remake=false, reload=true, + solver=FBDF(nlsolve=OrdinaryDiffEqNonlinearSolve.NLNewton(relax=0.8, max_iter=1000)) +) +OrdinaryDiffEqCore.set_proposed_dt!(s.integrator, dt) +sys = s.sys + +@info "System initialized at:" +toc() + +sys = s.sys +# # Stabilize system +s.integrator.ps[sys.stabilize] = true +for i in 1:10÷dt + next_step!(s; dt, vsm_interval=1) +end +s.integrator.ps[sys.stabilize] = false +plot(s, 0.0; zoom=true, front=false) + +# Function to step simulation with input u +function step_with_input_integ(x, u, _, p) + (s, set_x, set_ix, set_sx, sx, set_u, get_x, dt) = p + set_x(s.integrator, x) + set_sx(s.integrator, sx) + set_u(s.integrator, u) + OrdinaryDiffEqCore.set_t!(s.integrator, 0.0) + OrdinaryDiffEqCore.reinit!(s.integrator, s.integrator.u; reinit_dae=false) + OrdinaryDiffEqCore.set_proposed_dt!(s.integrator, dt) + OrdinaryDiffEqCore.step!(s.integrator, dt) + return get_x(s.integrator) +end + +# Get initial state +x_vec = KiteModels.get_nonstiff_unknowns(s) +sx_vec = KiteModels.get_stiff_unknowns(s) +set_ix = setu(s.integrator, Initial.(x_vec)) # set_ix might not be needed anymore if prob is removed +set_x = setu(s.integrator, x_vec) +set_sx = setu(s.integrator, sx_vec) +set_u = setu(s.integrator, collect(sys.set_values)) +get_x = getu(s.integrator, x_vec) +get_sx = getu(s.integrator, sx_vec) +sx = get_sx(s.integrator) +x0 = get_x(s.integrator) + +function test_response(s, input_range, input_idx, step_fn, u0, x_idxs=nothing; steps=1) + # If no x_idxs specified, default to angular velocities + if isnothing(x_idxs) + x_idxs = (length(x0)-8):(length(x0)-6) # Assuming ω_b are the last 3 non-stiff states before stiff ones + output_size = 3 + else + output_size = length(x_idxs) + end + + output = zeros(output_size, length(input_range)) + total_time = 0.0 + iter = 0 + + for (i, input_val) in enumerate(input_range) + u = copy(u0) + u[input_idx] += input_val + x = copy(x0) + set_ix(s.integrator, x) + OrdinaryDiffEqCore.reinit!(s.integrator) + sx_ = get_sx(s.integrator) + for _ in 1:steps # Use _ if i is not used inside the loop + p = (s, set_x, set_ix, set_sx, sx_, set_u, get_x, dt) # Pass set_ix even if unused by integ function + total_time += @elapsed x = step_fn(x, u, nothing, p) + iter += 1 + end + output[:, i] = x[x_idxs] + end + + times_rt = dt*iter/total_time + @info "Number of steps: $iter, Times realtime: $times_rt, Total time: $total_time" + return input_range, output, times_rt +end + +# Add helper function to find state indices +function find_state_index(x_vec, symbol) + # Compare the variables using isequal for symbolic equality + idx = findfirst(x -> isequal(x, symbol), x_vec) + isnothing(idx) && error("Symbol $symbol not found in state vector") + return idx +end + +function plot_input_output_relations(step_fn, suffix) + # Find relevant state indices + ω_idxs = [find_state_index(x_vec, sys.ω_b[i]) for i in 1:3] + twist_idx = find_state_index(sx_vec, sys.free_twist_angle[1]) + + # Test ranges + steer_range = range(-1, 1, length=100) + twist_range = range(-1, 1, length=100) + + # Test steering input vs omega + @info "Testing steering input response for $suffix..." + _, ω_steer_left, _ = test_response(s, steer_range, 2, step_fn, measure.set_values, ω_idxs) + _, ω_steer_right, _ = test_response(s, steer_range, 3, step_fn, measure.set_values, ω_idxs) + + # Test twist angle vs omega + @info "Testing twist angle response for $suffix..." + function step_with_twist(x, twist_val, _, p) + p[5][twist_idx] = twist_val[1] # Set twist angle directly in stiff state vector copy + return step_fn(x, measure.set_values, nothing, p) + end + _, ω_twist, _ = test_response(s, twist_range, 1, step_with_twist, zeros(3), ω_idxs) # u0 is zeros(3) here, seems ok + + return ω_steer_left, ω_steer_right, ω_twist, steer_range, twist_range +end + +ω_steer_left_integ, ω_steer_right_integ, ω_twist_integ, steer_range, twist_range = + plot_input_output_relations(step_with_input_integ, "integ") + +# --- Plot combined results --- +# Simplified data collection +steering_data = [ω_steer_left_integ[1,:], ω_steer_right_integ[1,:]] +steering_labels = ["Left Steering (integ)", "Right Steering (integ)"] +twist_data = [ω_twist_integ[1,:]] +twist_labels = ["Twist Input (integ)"] + +# Steering Plot +steering_plot = plotx(steer_range, + steering_data, # ω_b[1] data + [ω_steer_left_integ[2,:], ω_steer_right_integ[2,:]], # ω_b[2] data + [ω_steer_left_integ[3,:], ω_steer_right_integ[3,:]]; # ω_b[3] data + ylabels=["ω_b[1]", "ω_b[2]", "ω_b[3]"], + labels=[steering_labels, steering_labels, steering_labels], # Repeat labels for each subplot + fig="Steering Input vs Angular Velocity (Integrator)", + xlabel="Steering Input [Nm]") + +# Twist Plot +twist_plot = plotx(rad2deg.(twist_range), + twist_data, # ω_b[1] data + [ω_twist_integ[2,:]], # ω_b[2] data + [ω_twist_integ[3,:]]; # ω_b[3] data + ylabels=["ω_b[1]", "ω_b[2]", "ω_b[3]"], + labels=[twist_labels, twist_labels, twist_labels], # Repeat labels for each subplot + fig="Twist Angle vs Angular Velocity (Integrator)", + xlabel="Twist Angle [deg]") + +display(steering_plot) +display(twist_plot) diff --git a/examples/lin_ram_model.jl b/examples/lin_ram_model.jl new file mode 100644 index 000000000..e5bb3891a --- /dev/null +++ b/examples/lin_ram_model.jl @@ -0,0 +1,301 @@ +# Copyright (c) 2025 Bart van de Lint +# SPDX-License-Identifier: MPL-2.0 + +#= +This example demonstrates linearized model accuracy by comparing: +1. Nonlinear RamAirKite model simulation +2. Linearized state-space model simulation + +Both models start from the same operating point and are subjected +to identical steering inputs. The resulting state trajectories are +plotted together to visualize how well the linearized model +approximates the nonlinear dynamics. +=# + +using Timers +tic() +@info "Loading packages " + +using KiteModels, LinearAlgebra, Statistics, OrdinaryDiffEqCore +using ModelingToolkit + +PLOT = true +if PLOT + using Pkg + if ! ("LaTeXStrings" ∈ keys(Pkg.project().dependencies)) + using TestEnv; TestEnv.activate() + end + using ControlPlots, LaTeXStrings + import ControlPlots: plot +end +toc() + +include(joinpath(@__DIR__, "plotting.jl")) + +# TODO: use sparse autodiff + +# Simulation parameters +dt = 0.001 +total_time = 1.0 # Increased from 0.1s to 1.0s for better dynamics observation +vsm_interval = 3 +steps = Int(round(total_time / dt)) + +# Steering parameters +steering_freq = 1/2 # Hz - full left-right cycle frequency +steering_magnitude = 5.0 # Magnitude of steering input [Nm] + +# Initialize model +set = load_settings("system_ram.yaml") +set.segments = 3 +set_values = [-50.0, 0.0, 0.0] # Set values of the torques of the three winches. [Nm] +set.quasi_static = false +set.physical_model = "simple_ram" + +@info "Creating RamAirKite model..." +s = RamAirKite(set) +s.set.abs_tol = 1e-2 +s.set.rel_tol = 1e-2 +toc() + +measure = Measurement() + +# Define outputs for linearization - angular velocities +@variables ω_b(ModelingToolkit.t_nounits)[1:3] + +# Initialize at elevation with linearization outputs +s.point_system.winches[2].tether_length += 0.2 +s.point_system.winches[3].tether_length += 0.2 +measure.sphere_pos .= deg2rad.([65.0 65.0; 1.0 -1.0]) +KiteModels.init_sim!(s, measure; + remake=false, + reload=true, + lin_outputs=[ω_b...] # Specify which outputs to track in linear model +) +sys = s.sys + +@show rad2deg(s.integrator[sys.elevation]) + + +@info "System initialized at:" +toc() + +# --- Stabilize system at operating point --- +@info "Stabilizing system at operating point..." +s.integrator.ps[sys.stabilize] = true +stabilization_steps = Int(10 ÷ dt) +for i in 1:stabilization_steps + next_step!(s; dt, vsm_interval=0.05÷dt) +end +s.integrator.ps[sys.stabilize] = false + +# --- Linearize at operating point --- +@info "Linearizing system at operating point..." +@time (; A, B, C, D) = KiteModels.linearize(s) +@time (; A, B, C, D) = KiteModels.linearize(s) +@show norm(A) +@info "System linearized with matrix dimensions:" A=size(A) B=size(B) C=size(C) D=size(D) + +@show rad2deg(s.lin_prob[sys.elevation]) + +# --- Get operating point values --- +# Extract state and input at operating point +u_op = copy(s.integrator[sys.set_values]) +# Create a SysState to capture state values at operating point +sys_state_op = KiteModels.SysState(s) +# Also get the direct state vector for linear model +x_op = copy(s.integrator.u) + +# --- Create discretized system matrices --- +# Function to compute discretized system matrices using matrix exponential method +function discretize_linear_system(A, B, dt) + n = size(A, 1) + m = size(B, 2) + + # Create augmented matrix for simultaneous computation + M = [A B; zeros(m, n) zeros(m, m)] + + # Compute matrix exponential + expM = exp(M * dt) + + # Extract discretized matrices + Ad = expM[1:n, 1:n] + Bd = expM[1:n, n+1:n+m] + + return Ad, Bd +end + +# Discretize the continuous-time system for more accurate integration +@info "Discretizing system matrices..." +Ad, Bd = discretize_linear_system(A, B, dt) + +# Verify operating point stability +derivatives_at_op = A * zeros(size(A, 1)) + B * zeros(size(B, 2)) +@info "Derivatives at operating point (should be near zero):" norm(derivatives_at_op) + +# Create loggers +logger_nonlinear = Logger(length(s.point_system.points), steps) +logger_linear = deepcopy(logger_nonlinear) # Same structure for comparison + +# Initialize system state trackers +sys_state_nonlinear = KiteModels.SysState(s) +# For the linear model, we'll use a deviation from operating point +sys_state_linear = deepcopy(sys_state_op) + +# --- Prepare the simulation --- +simulation_time_points = zeros(Float64, steps) +# Pre-allocate arrays with fixed size instead of growing them dynamically +input_history = Vector{Vector{Float64}}(undef, steps) +perturbation_history = Vector{Vector{Float64}}(undef, steps) + +# --- Set up linear state tracking --- +# Create a linear model state vector that starts at zero (deviations from op point) +x_linear = zeros(length(x_op)) +# Create output vector for linear model +y_linear = zeros(size(C, 1)) + +@info "Starting side-by-side simulation..." +# Begin simulation with fixed number of steps +try + for i in 1:steps + global x_linear, sim_time + # Calculate current simulation time + sim_time = (i-1) * dt + simulation_time_points[i] = sim_time + + # --- Calculate time-varying steering inputs --- + # Use sinusoidal input for more realistic testing + steering = steering_magnitude * sin(2π * steering_freq * sim_time) + + # --- Nonlinear system simulation --- + # Compute control inputs: base winch force + steering + set_values_nonlinear = copy(u_op) + set_values_nonlinear .+= [0.0, steering, -steering] + input_history[i] = copy(set_values_nonlinear) + + # Compute perturbation from operating point + perturbation = set_values_nonlinear - u_op + perturbation_history[i] = copy(perturbation) + + # Step nonlinear simulation + (t_new, _) = next_step!(s, set_values_nonlinear; dt, vsm_interval=vsm_interval) + + # Log nonlinear state + KiteModels.update_sys_state!(sys_state_nonlinear, s) + log!(logger_nonlinear, sys_state_nonlinear) + + # --- Linear system simulation --- + # Use discretized state-space model + # x[k+1] = Ad*x[k] + Bd*u[k] + x_linear = Ad * x_linear + Bd * perturbation + + # Compute linear system output + y_linear .= C * x_linear + D * perturbation + sys_state_linear.turn_rates = sys_state_op.turn_rates .+ y_linear + + # Log linear state + log!(logger_linear, sys_state_linear) + end +catch e + if isa(e, AssertionError) + @show i + println(e) + else + rethrow(e) + end +end + +@info "Simulation completed:" +toc() + +# --- Save logs --- +save_log(logger_nonlinear, "nonlinear_model") +save_log(logger_linear, "linear_model") +lg_nonlinear = load_log("nonlinear_model") +lg_linear = load_log("linear_model") +sl_nonlinear = lg_nonlinear.syslog +sl_linear = lg_linear.syslog + +# --- Plot comparison results --- +# Extract necessary data +turn_rates_nonlinear_deg = rad2deg.(hcat(sl_nonlinear.turn_rates...)) +turn_rates_linear_deg = rad2deg.(hcat(sl_linear.turn_rates...)) + +# Format input data for plotting +steering_inputs = [inputs[2] - u_op[2] for inputs in input_history] + +# Create comparison plots using plotx +t_plot = sl_nonlinear.time + +# Prepare the data in the format expected by plotx +ω_x_comparison = [turn_rates_nonlinear_deg[1,:], turn_rates_linear_deg[1,:]] +ω_y_comparison = [turn_rates_nonlinear_deg[2,:], turn_rates_linear_deg[2,:]] +ω_z_comparison = [turn_rates_nonlinear_deg[3,:], turn_rates_linear_deg[3,:]] +steering_series = [steering_inputs] + +# Create the comparison plot +p_comparison = plotx(t_plot, + ω_x_comparison, + ω_y_comparison, + ω_z_comparison, + steering_series; + ylabels=["ω_x [°/s]", "ω_y [°/s]", "ω_z [°/s]", "Steering [Nm]"], + labels=[ + ["Nonlinear", "Linear"], + ["Nonlinear", "Linear"], + ["Nonlinear", "Linear"], + ["Input"] + ], + fig="Linear vs Nonlinear Model Comparison") +display(p_comparison) + +# --- Calculate error metrics --- +# Function to calculate normalized RMSE +function calculate_nrmse(actual, predicted) + # Trim to same length if needed + min_length = min(length(actual), length(predicted)) + # Calculate RMSE + rmse = sqrt(sum((actual[1:min_length] - predicted[1:min_length]).^2) / min_length) + # Normalize by range of actual values + range_actual = maximum(actual[1:min_length]) - minimum(actual[1:min_length]) + # Avoid division by zero + if abs(range_actual) < 1e-10 + return rmse # Return unnormalized if range is too small + else + return rmse / range_actual + end +end + +# Calculate both average L2 norm and NRMSE +len = turn_rates_nonlinear_deg[1,:] +error_ω_x = norm(turn_rates_nonlinear_deg[1,:] - turn_rates_linear_deg[1,:]) / len +error_ω_y = norm(turn_rates_nonlinear_deg[2,:] - turn_rates_linear_deg[2,:]) / len +error_ω_z = norm(turn_rates_nonlinear_deg[3,:] - turn_rates_linear_deg[3,:]) / len + +nrmse_ω_x = calculate_nrmse(turn_rates_nonlinear_deg[1,:], turn_rates_linear_deg[1,:]) +nrmse_ω_y = calculate_nrmse(turn_rates_nonlinear_deg[2,:], turn_rates_linear_deg[2,:]) +nrmse_ω_z = calculate_nrmse(turn_rates_nonlinear_deg[3,:], turn_rates_linear_deg[3,:]) + +@info "Error metrics (average L2 norm):" error_ω_x error_ω_y error_ω_z +@info "Normalized RMSE (lower is better):" nrmse_ω_x nrmse_ω_y nrmse_ω_z + +# --- Plot error over time --- +# Calculate difference between linear and nonlinear models +diff_ω_x = turn_rates_linear_deg[1,:] - turn_rates_nonlinear_deg[1,:] +diff_ω_y = turn_rates_linear_deg[2,:] - turn_rates_nonlinear_deg[2,:] +diff_ω_z = turn_rates_linear_deg[3,:] - turn_rates_nonlinear_deg[3,:] + +# Plot the differences +p_error = plotx(t_plot, + [diff_ω_x], + [diff_ω_y], + [diff_ω_z], + [steering_inputs]; + ylabels=["ω_x error [°/s]", "ω_y error [°/s]", "ω_z error [°/s]", "Steering [Nm]"], + labels=[ + ["Error"], + ["Error"], + ["Error"], + ["Input"] + ], + fig="Linear vs Nonlinear Model Error") +display(p_error) \ No newline at end of file diff --git a/examples/plotting.jl b/examples/plotting.jl index fb702854d..0b1d356f2 100644 --- a/examples/plotting.jl +++ b/examples/plotting.jl @@ -1,16 +1,16 @@ # Copyright (c) 2022, 2024 Uwe Fechner # SPDX-License-Identifier: MIT -function plot(sys::PointMassSystem, reltime; kite_pos=nothing, zoom=false, front=false) +function plot(sys::PointMassSystem, reltime; kite_pos=nothing, e_z=zeros(3), zoom=false, front=false) pos = [sys.points[i].pos_w for i in eachindex(sys.points)] !isnothing(kite_pos) && (pos = [pos..., kite_pos]) seg = [[sys.segments[i].points[1], sys.segments[i].points[2]] for i in eachindex(sys.segments)] if zoom && !front - xlim = (pos[end][1] - 5, pos[end][1]+5) - ylim = (pos[end][3] - 8, pos[end][3]+2) + xlim = (pos[end][1] - 6, pos[end][1]+6) + ylim = (pos[end][3] - 10, pos[end][3]+2) elseif zoom && front - xlim = (pos[end][2] - 5, pos[end][2]+5) - ylim = (pos[end][3] - 8, pos[end][3]+2) + xlim = (pos[end][2] - 6, pos[end][2]+6) + ylim = (pos[end][3] - 10, pos[end][3]+2) elseif !zoom && !front xlim = (0, 60) ylim = (0, 60) @@ -27,5 +27,6 @@ function plot(s::RamAirKite, reltime; kwargs...) for (i, point) in enumerate(s.point_system.points) point.pos_w .= pos[:, i] end - plot(s.point_system, reltime; kite_pos, kwargs...) + e_z = s.integrator[s.sys.e_z] + plot(s.point_system, reltime; kite_pos, e_z, kwargs...) end \ No newline at end of file diff --git a/examples/ram_air_kite.jl b/examples/ram_air_kite.jl index 9f99bd8a3..a7aa71ec7 100644 --- a/examples/ram_air_kite.jl +++ b/examples/ram_air_kite.jl @@ -5,7 +5,7 @@ using Timers tic() @info "Loading packages " -using KiteModels, LinearAlgebra +using KiteModels, LinearAlgebra, Statistics PLOT = true if PLOT @@ -24,49 +24,55 @@ include(joinpath(@__DIR__, "plotting.jl")) # Simulation parameters dt = 0.05 total_time = 10 # Longer simulation to see oscillations -vsm_interval = 2 +vsm_interval = 3 steps = Int(round(total_time / dt)) # Steering parameters steering_freq = 1/2 # Hz - full left-right cycle frequency -steering_magnitude = 5.0 # Magnitude of steering input [Nm] +steering_magnitude = 10.0 # Magnitude of steering input [Nm] # Initialize model -set = se("system_ram.yaml") +set = load_settings("system_ram.yaml") set.segments = 3 set_values = [-50, 0.0, 0.0] # Set values of the torques of the three winches. [Nm] set.quasi_static = false +set.physical_model = "ram" @info "Creating wing, aero, vsm_solver, point_system and s:" -wing = RamAirWing(set; prn=false) -aero = BodyAerodynamics([wing]) -vsm_solver = Solver(aero; solver_type=NONLIN, atol=2e-8, rtol=2e-8) -point_system = PointMassSystem(set, wing) -s = RamAirKite(set, aero, vsm_solver, point_system) +s = RamAirKite(set) +s.set.abs_tol = 1e-2 +s.set.rel_tol = 1e-2 toc() +# init_Q_b_w, R_b_w = KiteModels.measure_to_q(measure) +# init_kite_pos = init!(s.point_system, s.set, R_b_w, init_Q_b_w) +# plot(s.point_system, 0.0; zoom=false, front=true) + measure = Measurement() -s.set.abs_tol = 1e-5 -s.set.rel_tol = 1e-4 # Initialize at elevation -measure.sphere_pos .= deg2rad.([60.0 60.0; 1.0 -1.0]) -KiteModels.init_sim!(s, measure) +s.point_system.winches[2].tether_length += 0.2 +s.point_system.winches[3].tether_length += 0.2 +measure.sphere_pos .= deg2rad.([70.0 70.0; 1.0 -1.0]) +KiteModels.init_sim!(s, measure; remake=false, reload=true) sys = s.sys @info "System initialized at:" toc() -# Stabilize system -s.integrator.ps[sys.steady] = true -next_step!(s; dt=10.0, vsm_interval=1) -s.integrator.ps[sys.steady] = false +# # Stabilize system +s.integrator.ps[sys.stabilize] = true +for i in 1:10÷dt + next_step!(s; dt, vsm_interval=1) +end +s.integrator.ps[sys.stabilize] = false logger = Logger(length(s.point_system.points), steps) sys_state = KiteModels.SysState(s) t = 0.0 runtime = 0.0 integ_runtime = 0.0 +bias = 2.0 try while t < total_time @@ -77,43 +83,27 @@ try # Calculate steering inputs based on cosine wave steering = steering_magnitude * cos(2π * steering_freq * t+0.1) set_values = -s.set.drum_radius .* s.integrator[sys.winch_force] + _vsm_interval = 1 if t > 1.0 set_values .+= [0.0, steering, -steering] # Opposite steering for left/right + _vsm_interval = vsm_interval + else + set_values[2] += bias end - + # Step simulation - steptime = @elapsed (t_new, integ_steptime) = next_step!(s, set_values; dt, vsm_interval) + steptime = @elapsed (t_new, integ_steptime) = next_step!(s, set_values; dt, vsm_interval=_vsm_interval) t = t_new - 10.0 # Adjust for initial stabilization time - + # Track performance after initial transient if (t > total_time/2) runtime += steptime integ_runtime += integ_steptime + s.integrator.ps[sys.twist_damp] = 10 end - + # Log state variables KiteModels.update_sys_state!(sys_state, s) - sys_state.var_01 = s.integrator[sys.ω_b[1]] - sys_state.var_02 = s.integrator[sys.ω_b[2]] - sys_state.var_03 = s.integrator[sys.ω_b[3]] - - sys_state.var_04 = s.integrator[sys.tether_vel[2]] - sys_state.var_05 = s.integrator[sys.tether_vel[3]] - - sys_state.var_06 = s.integrator[sys.aero_force_b[3]] - sys_state.var_07 = s.integrator[sys.aero_moment_b[2]] - sys_state.var_08 = s.integrator[sys.group_aero_moment[1]] - - sys_state.var_09 = s.integrator[sys.twist_angle[1]] - sys_state.var_10 = s.integrator[sys.twist_angle[2]] - sys_state.var_11 = s.integrator[sys.twist_angle[3]] - sys_state.var_12 = s.integrator[sys.twist_angle[4]] - - sys_state.var_13 = s.integrator[sys.pulley_l0[1]] - sys_state.var_14 = s.integrator[sys.pulley_l0[2]] - - sys_state.var_15 = rad2deg(calc_aoa(s)) - log!(logger, sys_state) end catch e @@ -133,23 +123,31 @@ save_log(logger, "tmp") lg =load_log("tmp") sl = lg.syslog -p = plotx(sl.time .- 10, - [rad2deg.(sl.var_01), rad2deg.(sl.var_02), rad2deg.(sl.var_03)], - [c(sl.var_04), c(sl.var_05)], - [c(sl.var_06), c(sl.var_07), c(sl.var_08)], - [rad2deg.(c(sl.var_09)), rad2deg.(c(sl.var_10)), rad2deg.(c(sl.var_11)), rad2deg.(c(sl.var_12))], - [c(sl.var_13), c(sl.var_14)], - [c(sl.var_15)], - [rad2deg.(c(sl.heading))]; - ylabels=["turn rates [°/s]", L"v_{ro}~[m/s]", "vsm", "twist [°]", "pulley", "AoA [°]", "heading [°]"], +# --- Updated Plotting --- +# Extract necessary data using meaningful names +turn_rates_deg = rad2deg.(hcat(sl.turn_rates...)) +v_reelout_23 = [sl.v_reelout[i][2] for i in eachindex(sl.v_reelout)], [sl.v_reelout[i][3] for i in eachindex(sl.v_reelout)] # Winch 2 and 3 +aero_force_z = [sl.aero_force_b[i][3] for i in eachindex(sl.aero_force_b)] +aero_moment_y = [sl.aero_moment_b[i][2] for i in eachindex(sl.aero_moment_b)] +twist_angles_deg = rad2deg.(hcat(sl.twist_angles...)) +AoA_deg = rad2deg.(sl.AoA) +heading_deg = rad2deg.(sl.heading) + +p = plotx(sl.time .- 10, + [turn_rates_deg[1,:], turn_rates_deg[2,:], turn_rates_deg[3,:]], + v_reelout_23, + [aero_force_z, aero_moment_y], + [twist_angles_deg[1,:], twist_angles_deg[2,:], twist_angles_deg[3,:], twist_angles_deg[4,:]], + [AoA_deg], + [heading_deg]; + ylabels=["turn rates [°/s]", L"v_{ro}~[m/s]", "aero F/M", "twist [°]", "AoA [°]", "heading [°]"], ysize=10, labels=[ - [L"ω_x", L"ω_y", L"ω_z"], - ["vel[1]", "vel[2]"], - ["force[3]", "kite moment[2]", "group moment[1]"], - ["twist_angle[1]", "twist_angle[2]", "twist_angle[3]", "twist_angle[4]"], - ["pulley_l0[1]", "pulley_l0[2]"], - ["angle of attack"], + [L"\omega_x", L"\omega_y", L"\omega_z"], + ["v_ro[2]", "v_ro[3]"], + [L"F_{aero,z}", L"M_{aero,y}"], + ["twist[1]", "twist[2]", "twist[3]", "twist[4]"], + ["AoA"], ["heading"] ], fig="Oscillating Steering Input Response") diff --git a/examples/simple_ram_air_kite.jl b/examples/simple_ram_air_kite.jl new file mode 100644 index 000000000..ddfab861e --- /dev/null +++ b/examples/simple_ram_air_kite.jl @@ -0,0 +1,140 @@ +# Copyright (c) 2025 Bart van de Lint +# SPDX-License-Identifier: MPL-2.0 + +using KiteModels, LinearAlgebra +using Pkg +if ! ("LaTeXStrings" ∈ keys(Pkg.project().dependencies)) + using TestEnv; TestEnv.activate() +end +using LaTeXStrings + +PLOT = true +if PLOT + using ControlPlots +end + +include(joinpath(@__DIR__, "plotting.jl")) + +# Simulation parameters +dt = 0.05 +total_time = 10 # Longer simulation to see oscillations +vsm_interval = 2 +steps = Int(round(total_time / dt)) + +# Steering parameters +steering_freq = 1/2 # Hz - full left-right cycle frequency +steering_magnitude = 1.0 # Magnitude of steering input [Nm] + +# Initialize model +set = se("system_ram.yaml") +set.segments = 2 +set.quasi_static = true +set.bridle_fracs = [0.0, 0.93] +set.physical_model = "simple_ram" + +wing = RamAirWing(set; prn=false, n_groups=2) +aero = BodyAerodynamics([wing]) +vsm_solver = Solver(aero; solver_type=NONLIN, atol=2e-8, rtol=2e-8) +point_system = create_simple_ram_point_system(set, wing) +s = RamAirKite(set, aero, vsm_solver, point_system) + +measure = Measurement() +measure.set_values .= [-55, -4.0, -4.0] # Set values of the torques of the three winches. [Nm] +set_values = measure.set_values +s.set.abs_tol = 1e-5 +s.set.rel_tol = 1e-4 + +# Initialize at elevation +measure.sphere_pos .= deg2rad.([83.0 83.0; 1.0 -1.0]) +KiteModels.init_sim!(s, measure; remake=false) +sys = s.sys + +# Stabilize system +s.integrator.ps[sys.steady] = true +next_step!(s; dt=10.0, vsm_interval=1) +s.integrator.ps[sys.steady] = false + +logger = Logger(length(s.point_system.points), steps) +sys_state = KiteModels.SysState(s) +t = 0.0 +runtime = 0.0 +integ_runtime = 0.0 + +try + while t < total_time + global t, set_values, runtime, integ_runtime + PLOT && plot(s, t; zoom=false, front=false) + + # Calculate steering inputs based on cosine wave + steering = steering_magnitude * cos(2π * steering_freq * t) + set_values = -s.set.drum_radius .* s.integrator[sys.winch_force] + set_values .+= [0.0, steering, -steering] # Opposite steering for left/right + + # Step simulation + steptime = @elapsed (t_new, integ_steptime) = next_step!(s, set_values; dt, vsm_interval) + t = t_new - 10.0 # Adjust for initial stabilization time + + # Track performance after initial transient + if (t > total_time/2) + runtime += steptime + integ_runtime += integ_steptime + end + + # Log state variables + KiteModels.update_sys_state!(sys_state, s) + sys_state.var_01 = s.integrator[sys.ω_b[1]] + sys_state.var_02 = s.integrator[sys.ω_b[2]] + sys_state.var_03 = s.integrator[sys.ω_b[3]] + + sys_state.var_04 = s.integrator[sys.tether_vel[2]] + sys_state.var_05 = s.integrator[sys.tether_vel[3]] + + sys_state.var_06 = s.integrator[sys.aero_force_b[3]] + sys_state.var_07 = s.integrator[sys.aero_moment_b[2]] + sys_state.var_08 = s.integrator[sys.group_aero_moment[1]] + + sys_state.var_09 = s.integrator[sys.twist_angle[1]] + sys_state.var_10 = s.integrator[sys.twist_angle[2]] + + sys_state.var_11 = s.integrator[sys.angle_of_attack] # TODO: investigate why different from vsm aoa + + log!(logger, sys_state) + end +catch e + if isa(e, AssertionError) + @show t + println(e) + else + rethrow(e) + end +end + +# Plot results +c = collect +save_log(logger, "tmp") +lg =load_log("tmp") +sl = lg.syslog + +if PLOT + p = plotx(sl.time .- 10, + [rad2deg.(sl.var_01), rad2deg.(sl.var_02), rad2deg.(sl.var_03)], + [c(sl.var_04), c(sl.var_05)], + [c(sl.var_06), c(sl.var_07), c(sl.var_08)], + [rad2deg.(c(sl.var_09)), rad2deg.(c(sl.var_10))], + [rad2deg.(c(sl.var_11))], + [rad2deg.(c(sl.heading))]; + ylabels=["turn rates [°/s]", L"v_{ro}~[m/s]", "vsm", "twist [°]", "AoA [°]", "heading [°]"], + ysize=10, + labels=[ + [L"ω_x", L"ω_y", L"ω_z"], + ["vel[1]", "vel[2]"], + ["force[3]", "kite moment[2]", "group moment[1]"], + ["twist_angle[1]", "twist_angle[2]"], + ["angle of attack"], + ["heading"] + ], + fig="Oscillating Steering Input Response") + display(p) +end + +@info "Performance:" times_realtime=(total_time/2)/runtime integrator_times_realtime=(total_time/2)/integ_runtime diff --git a/examples/test_init_1p.jl b/examples/test_init_1p.jl index 8046edb21..8e2bb0351 100644 --- a/examples/test_init_1p.jl +++ b/examples/test_init_1p.jl @@ -34,7 +34,7 @@ DEPOWER = 0.47:-0.005:0.355 elev = set.elevation i = 1 set.v_wind = V_WIND # 25 -logger::Logger = Logger(set.segments + 1, STEPS) +logger = Logger(set.segments + 1, STEPS) kcu::KCU = KCU(set) kps3::KPS3 = KPS3(kcu) diff --git a/examples/test_init_4p.jl b/examples/test_init_4p.jl index 0f894f22d..7df06304d 100644 --- a/examples/test_init_4p.jl +++ b/examples/test_init_4p.jl @@ -35,7 +35,7 @@ UPWIND_DIR = -pi/2 +deg2rad(10) elev = set.elevation i = 1 set.v_wind = V_WIND # 25 -logger::Logger = Logger(set.segments + 5, STEPS) +logger = Logger(set.segments + 5, STEPS) kcu::KCU = KCU(set) kps4::KPS4 = KPS4(kcu) diff --git a/src/KiteModels.jl b/src/KiteModels.jl index af173d8d9..d7ade9e26 100644 --- a/src/KiteModels.jl +++ b/src/KiteModels.jl @@ -45,6 +45,7 @@ export calc_azimuth_north, calc_azimuth_east export winch_force, lift_drag, cl_cd, lift_over_drag, unstretched_length, tether_length, v_wind_kite # getters export calculate_rotational_inertia! export kite_ref_frame, orient_euler, spring_forces, upwind_dir, copy_model_settings, menu2 +export create_ram_point_system, create_simple_ram_point_system import LinearAlgebra: norm set_zero_subnormals(true) # required to avoid drastic slow down on Intel CPUs when numbers become very small @@ -69,6 +70,7 @@ const SimFloat = Float64 Basic 3-dimensional vector, stack allocated, mutable. """ const KVec3 = MVector{3, SimFloat} +const KVec4 = MVector{4, SimFloat} """ const SVec3 = SVector{3, SimFloat} @@ -99,11 +101,11 @@ function __init__() end include("KPS4.jl") # include code, specific for the four point kite model -include("ram_air_kite.jl") # include code, specific for the four point 3 line kite model +include("point_mass_system.jl") +include("ram_air_kite.jl") # include code, specific for the ram air kite model include("mtk_model.jl") include("KPS3.jl") # include code, specific for the one point kite model include("init.jl") # functions to calculate the initial state vector, the initial masses and initial springs -include("point_mass_system.jl") function menu2() Main.include("examples/menu2.jl") diff --git a/src/mtk_model.jl b/src/mtk_model.jl index a4a086e84..b4689d344 100644 --- a/src/mtk_model.jl +++ b/src/mtk_model.jl @@ -1,5 +1,5 @@ -# Copyright (c) 2024 Bart van de Lint -# SPDX-License-Identifier: MIT +# Copyright (c) 2025 Bart van de Lint +# SPDX-License-Identifier: MPL-2.0 # ==================== mtk model functions ================================================ # Implementation of the ram air kite model using ModelingToolkit.jl @@ -10,36 +10,17 @@ end function calc_moment_acc(winch::TorqueControlledMachine, tether_vel, norm_, set_torque) calc_acceleration(winch, tether_vel, norm_; set_speed=nothing, set_torque, use_brake=false) end -@register_symbolic calc_speed_acc(winch::AsyncMachine, tether_vel, norm_, set_speed) -@register_symbolic calc_moment_acc(winch::TorqueControlledMachine, tether_vel, norm_, set_torque) -function sym_interp(interp::Function, aoa, trailing_edge_angle) - return interp(rad2deg(aoa), rad2deg(trailing_edge_angle-aoa)) # TODO: register callable struct https://docs.sciml.ai/Symbolics/dev/manual/functions/#Symbolics.@register_array_symbolic +function calc_heading(e_x) + return atan(e_x[2], e_x[1]) end -@register_symbolic sym_interp(interp::Function, aoa, trailing_edge_angle) -function sym_normalize(vec) - return vec / norm(vec) -end -@register_symbolic sym_normalize(vec) - -function rotate_by_quaternion(q, v) - p = [0, v...] - return quaternion_multiply(q, quaternion_multiply(p, quaternion_conjugate(q)))[2:4] +function calc_angle_of_attack(va_kite_b) + return atan(va_kite_b[3], va_kite_b[1]) end -function quaternion_conjugate(q) - return [q[1], -q[2], -q[3], -q[4]] -end - -function quaternion_multiply(q1, q2) - w1, x1, y1, z1 = q1 - w2, x2, y2, z2 = q2 - w = w1*w2 - x1*x2 - y1*y2 - z1*z2 - x = w1*x2 + x1*w2 + y1*z2 - z1*y2 - y = w1*y2 - x1*z2 + y1*w2 + z1*x2 - z = w1*z2 + x1*y2 - y1*x2 + z1*w2 - return [w, x, y, z] +function sym_normalize(vec) + return vec / norm(vec) end function quaternion_to_rotation_matrix(q) @@ -102,7 +83,7 @@ pulley dynamics and winch forces. - `wind_vec_gnd`: Ground wind vector - `group_aero_moment`: Aerodynamic moments per group - `twist_angle`: Twist angles per group -- `steady`: Whether in steady mode +- `stabilize`: Whether in stabilize mode # Returns Tuple containing: @@ -113,7 +94,7 @@ Tuple containing: - Tether moments on kite """ function force_eqs!(s, system, eqs, defaults, guesses; - R_b_w, kite_pos, kite_vel, wind_vec_gnd, group_aero_moment, twist_angle, steady) + R_b_w, kite_pos, kite_vel, wind_vec_gnd, group_aero_moment, twist_angle, twist_ω, stabilize, set_values, fix_nonstiff) @parameters acc_multiplier = 1 @@ -187,23 +168,31 @@ function force_eqs!(s, system, eqs, defaults, guesses; found += 1 end end - !(found == 1) && throw(ArgumentError("Kite point number $(point.idx) is part of $found groups, - and should be part of exactly 1 groups.")) + !(found in [0,1]) && throw(ArgumentError("Kite point number $(point.idx) is part of $found groups, + and should be part of exactly 0 or 1 groups.")) - group = groups[group_idx] - if point.idx == group.points[group.fixed_index] - pos_b = point.pos_b - else - fixed_pos = points[group.points[group.fixed_index]].pos_b + if found == 1 + group = groups[group_idx] + fixed_pos = group.le_pos chord_b = point.pos_b - fixed_pos normal = chord_b × group.y_airf pos_b = fixed_pos + cos(twist_angle[group_idx]) * chord_b - sin(twist_angle[group_idx]) * normal + eqs = [ + eqs + tether_kite_moment[:, point.idx] ~ zeros(3) + ] + else + pos_b = point.pos_b + eqs = [ + eqs + tether_r[:, point.idx] ~ pos[:, point.idx] - kite_pos + tether_kite_moment[:, point.idx] ~ tether_r[:, point.idx] × tether_kite_force[:, point.idx] + ] end + eqs = [ eqs tether_kite_force[:, point.idx] ~ point_force[:, point.idx] - tether_r[:, point.idx] ~ pos[:, point.idx] - kite_pos - tether_kite_moment[:, point.idx] ~ tether_r[:, point.idx] × point_force[:, point.idx] pos[:, point.idx] ~ kite_pos + R_b_w * pos_b vel[:, point.idx] ~ zeros(3) acc[:, point.idx] ~ zeros(3) @@ -220,7 +209,7 @@ function force_eqs!(s, system, eqs, defaults, guesses; n = sym_normalize(kite_pos) n = n * (p ⋅ n) r = (p - n) # https://en.wikipedia.org/wiki/Distance_from_a_point_to_a_line#Vector_formulation - @parameters bridle_damp = 100 + @parameters bridle_damp = 1.0 @parameters measured_ω_z = 0.6 eqs = [ eqs @@ -228,7 +217,7 @@ function force_eqs!(s, system, eqs, defaults, guesses; D(vel[:, point.idx]) ~ acc_multiplier * acc[:, point.idx] - in_bridle * bridle_damp * (vel[:, point.idx] - kite_vel) acc[:, point.idx] ~ point_force[:, point.idx] / mass + [0, 0, -G_EARTH] + - ifelse.(steady==true, r * norm(measured_ω_z)^2, zeros(3)) # TODO: add other steady accelerations + ifelse.(stabilize==true, r * norm(measured_ω_z)^2, zeros(3)) # TODO: add other stabilize accelerations ] defaults = [ defaults @@ -259,10 +248,10 @@ function force_eqs!(s, system, eqs, defaults, guesses; trailing_edge_ω(t)[eachindex(groups)] # angular rate trailing_edge_α(t)[eachindex(groups)] # angular acc free_twist_angle(t)[eachindex(groups)] - twist_ω(t)[eachindex(groups)] # angular rate twist_α(t)[eachindex(groups)] # angular acc group_tether_moment(t)[eachindex(groups)] - tether_moment(t)[eachindex(groups), 1:length(groups[1].points)-1] + tether_force(t)[eachindex(groups), eachindex(groups[1].points)] + tether_moment(t)[eachindex(groups), eachindex(groups[1].points)] end for group in groups @@ -270,18 +259,19 @@ function force_eqs!(s, system, eqs, defaults, guesses; x_airf = normalize(group.chord) init_z_airf = x_airf × group.y_airf z_airf = x_airf * sin(twist_angle[group.idx]) + init_z_airf * cos(twist_angle[group.idx]) - moving_points = filter(p -> p != group.points[group.fixed_index], group.points) - for (i, point_idx) in enumerate(moving_points) - r = (points[point_idx].pos_b - points[group.points[group.fixed_index]].pos_b) ⋅ normalize(group.chord) + for (i, point_idx) in enumerate(group.points) + r = (points[point_idx].pos_b - (group.le_pos + group.moment_frac*group.chord)) ⋅ normalize(group.chord) eqs = [ eqs - tether_moment[group.idx, i] ~ r * (point_force[:, point_idx] ⋅ (R_b_w * -z_airf)) + tether_force[group.idx, i] ~ (point_force[:, point_idx] ⋅ (R_b_w * -z_airf)) + tether_moment[group.idx, i] ~ r * tether_force[group.idx, i] ] end inertia = 1/3 * (s.set.mass/length(groups)) * (norm(group.chord))^2 # plate inertia around leading edge @assert !(inertia ≈ 0.0) - @parameters twist_damp = s.set.quasi_static ? 200 : 100 + @parameters twist_damp = 50 + eqs = [ eqs group_tether_moment[group.idx] ~ sum(tether_moment[group.idx, :]) @@ -291,8 +281,8 @@ function force_eqs!(s, system, eqs, defaults, guesses; if group.type == DYNAMIC eqs = [ eqs - D(free_twist_angle[group.idx]) ~ twist_ω[group.idx] - D(twist_ω[group.idx]) ~ twist_α[group.idx] - twist_damp * twist_ω[group.idx] + D(free_twist_angle[group.idx]) ~ ifelse(fix_nonstiff==true, 0, twist_ω[group.idx]) + D(twist_ω[group.idx]) ~ ifelse(fix_nonstiff==true, 0, twist_α[group.idx] - twist_damp * twist_ω[group.idx]) ] defaults = [ defaults @@ -341,10 +331,12 @@ function force_eqs!(s, system, eqs, defaults, guesses; end for segment in segments p1, p2 = segment.points[1], segment.points[2] - guesses = [ - guesses - [segment_vec[i, segment.idx] => points[p2].pos_w[i] - points[p1].pos_w[i] for i in 1:3] - ] + if s.set.quasi_static + guesses = [ + guesses + [segment_vec[i, segment.idx] => points[p2].pos_w[i] - points[p1].pos_w[i] for i in 1:3] + ] + end if segment.type == BRIDLE in_pulley = 0 @@ -401,15 +393,18 @@ function force_eqs!(s, system, eqs, defaults, guesses; end stiffness_m = s.set.e_tether * (segment.diameter/2)^2 * pi - (segment.type == BRIDLE) && (compression_frac = 0.1) - (segment.type == POWER) && (compression_frac = 0.1) - (segment.type == STEERING) && (compression_frac = 0.1) - - @parameters stiffness_frac = 0.1 + @parameters stiffness_frac = 0.01 (segment.type == BRIDLE) && (stiffness_m = stiffness_frac * stiffness_m) damping_m = (s.set.damping / s.set.c_spring) * stiffness_m + if segment.compression_frac ≈ 1.0 + eqs = [eqs; stiffness[segment.idx] ~ stiffness_m / len[segment.idx]] + else + eqs = [eqs; stiffness[segment.idx] ~ ifelse(len[segment.idx] > l0[segment.idx], + stiffness_m / len[segment.idx], + segment.compression_frac * stiffness_m / len[segment.idx])] + end eqs = [ eqs # spring force equations @@ -418,10 +413,6 @@ function force_eqs!(s, system, eqs, defaults, guesses; unit_vec[:, segment.idx] ~ segment_vec[:, segment.idx]/len[segment.idx] rel_vel[:, segment.idx] ~ vel[:, p1] - vel[:, p2] spring_vel[segment.idx] ~ rel_vel[:, segment.idx] ⋅ unit_vec[:, segment.idx] - stiffness[segment.idx] ~ ifelse(len[segment.idx] > l0[segment.idx], - stiffness_m / len[segment.idx], - compression_frac * stiffness_m / len[segment.idx] - ) damping[segment.idx] ~ damping_m / len[segment.idx] spring_force[segment.idx] ~ (stiffness[segment.idx] * (len[segment.idx] - l0[segment.idx]) - damping[segment.idx] * spring_vel[segment.idx]) @@ -448,7 +439,7 @@ function force_eqs!(s, system, eqs, defaults, guesses; pulley_force(t)[eachindex(pulleys)] pulley_acc(t)[eachindex(pulleys)] end - @parameters pulley_damp = 20 + @parameters pulley_damp = 1.0 for pulley in pulleys segment = segments[pulley.segments[1]] mass_per_meter = s.set.rho_tether * π * (segment.diameter/2)^2 @@ -486,7 +477,6 @@ function force_eqs!(s, system, eqs, defaults, guesses; end # ==================== WINCHES ==================== # - @parameters set_values(t)[eachindex(winches)] = zeros(length(winches)) @variables begin tether_vel(t)[eachindex(winches)] tether_acc(t)[eachindex(winches)] @@ -500,14 +490,14 @@ function force_eqs!(s, system, eqs, defaults, guesses; end eqs = [ eqs - D(tether_length[winch.idx]) ~ tether_vel[winch.idx] - D(tether_vel[winch.idx]) ~ ifelse(steady==true, 0, tether_acc[winch.idx]) + D(tether_length[winch.idx]) ~ ifelse(fix_nonstiff==true, 0, ifelse(stabilize==true, 0, tether_vel[winch.idx])) + D(tether_vel[winch.idx]) ~ ifelse(fix_nonstiff==true, 0, ifelse(stabilize==true, 0, tether_acc[winch.idx])) tether_acc[winch.idx] ~ calc_moment_acc( # TODO: moment and speed control winch.model, tether_vel[winch.idx], winch_force[winch.idx], - set_values[winch.idx] - ) + set_values[winch.idx] # TODO: SISO mode to test linear model accuracy, MHE to check if it is better, use stabilize as state estimator + ) # TODO: add brakes and P. This makes SISO system easy winch_force[winch.idx] ~ F ] defaults = [ @@ -538,20 +528,20 @@ angular velocities and accelerations, and forces/moments. - `kite_pos`: Kite position vector - `kite_vel`: Kite velocity vector - `kite_acc`: Kite acceleration vector -- `init_Q_b_w`: Initial quaternion orientation -- `init_kite_pos`: Initial kite position -- `steady`: Whether in steady mode +- `stabilize`: Whether in stabilize mode # Returns Tuple of updated equations and defaults """ function diff_eqs!(s, eqs, defaults; tether_kite_force, tether_kite_moment, aero_force_b, - aero_moment_b, ω_b, R_b_w, kite_pos, kite_vel, kite_acc, init_Q_b_w, init_kite_pos, steady + aero_moment_b, ω_b, α_b, R_b_w, kite_pos, kite_vel, kite_acc, stabilize, fix_nonstiff ) + kite = s.point_system.kite @variables begin # potential differential variables kite_acc_b(t)[1:3] - α_b(t)[1:3] # angular acceleration in principal frame + α_b_damped(t)[1:3] + ω_b_stable(t)[1:3] # rotations and frames Q_b_w(t)[1:4] # quaternion orientation of the kite body frame relative to the world frame @@ -564,72 +554,95 @@ function diff_eqs!(s, eqs, defaults; tether_kite_force, tether_kite_moment, aero end @parameters ω_damp = 150 - Ω = [0 -ω_b[1] -ω_b[2] -ω_b[3]; - ω_b[1] 0 ω_b[3] -ω_b[2]; - ω_b[2] -ω_b[3] 0 ω_b[1]; - ω_b[3] ω_b[2] -ω_b[1] 0] + Ω(ω) = [0 -ω[1] -ω[2] -ω[3]; + ω[1] 0 ω[3] -ω[2]; + ω[2] -ω[3] 0 ω[1]; + ω[3] ω[2] -ω[1] 0] I_b = [s.wing.inertia_tensor[1,1], s.wing.inertia_tensor[2,2], s.wing.inertia_tensor[3,3]] @assert !(s.set.mass ≈ 0) + axis = sym_normalize(kite_pos) + axis_b = R_b_w' * axis eqs = [ eqs [D(Q_b_w[i]) ~ Q_vel[i] for i in 1:4] - [Q_vel[i] ~ 0.5 * sum(Ω[i, j] * Q_b_w[j] for j in 1:4) for i in 1:4] - [R_b_w[:, i] ~ quaternion_to_rotation_matrix(Q_b_w)[:, i] for i in 1:3] - - D(ω_b[1]) ~ α_b[1] - D(ω_b[2]) ~ α_b[2] - ω_damp * ω_b[2] - D(ω_b[3]) ~ ifelse(steady==true, 0, α_b[3]) + [Q_vel[i] ~ 0.5 * sum(Ω(ω_b_stable)[i, j] * Q_b_w[j] for j in 1:4) for i in 1:4] + ω_b_stable ~ ifelse.(fix_nonstiff==true, + zeros(3), + ifelse.(stabilize==true, + ω_b - ω_b ⋅ axis_b * axis_b, + ω_b + ) + ) + D(ω_b) ~ ifelse.(fix_nonstiff==true, + zeros(3), + ifelse.(stabilize==true, + α_b_damped - α_b_damped ⋅ axis_b * axis_b, + α_b_damped + ) + ) + α_b_damped ~ [α_b[1], α_b[2] - ω_damp*ω_b[2], α_b[3]] + [R_b_w[:, i] ~ quaternion_to_rotation_matrix(Q_b_w)[:, i] for i in 1:3] + α_b[1] ~ (moment_b[1] + (I_b[2] - I_b[3]) * ω_b[2] * ω_b[3]) / I_b[1] α_b[2] ~ (moment_b[2] + (I_b[3] - I_b[1]) * ω_b[3] * ω_b[1]) / I_b[2] α_b[3] ~ (moment_b[3] + (I_b[1] - I_b[2]) * ω_b[1] * ω_b[2]) / I_b[3] total_tether_kite_moment ~ [sum(tether_kite_moment[i, :]) for i in 1:3] moment_b ~ aero_moment_b + R_b_w' * total_tether_kite_moment - D(kite_pos) ~ kite_vel - D(kite_vel) ~ kite_acc - kite_acc ~ R_b_w * [ifelse(steady==true, 0, kite_acc_b[1]), - ifelse(steady==true, 0, kite_acc_b[2]), - kite_acc_b[3]] - kite_acc_b ~ (R_b_w' * total_tether_kite_force + aero_force_b) / s.set.mass + D(kite_pos) ~ ifelse.(fix_nonstiff==true, + zeros(3), + ifelse.(stabilize==true, + kite_vel ⋅ axis * axis, + kite_vel + ) + ) + D(kite_vel) ~ ifelse.(fix_nonstiff==true, + zeros(3), + ifelse.(stabilize==true, + kite_acc ⋅ axis * axis, + kite_acc + ) + ) + kite_acc ~ (total_tether_kite_force + R_b_w * aero_force_b) / s.set.mass total_tether_kite_force ~ [sum(tether_kite_force[i, :]) for i in 1:3] ] defaults = [ defaults - [Q_b_w[i] => init_Q_b_w[i] for i in 1:4] - [ω_b[i] => 0 for i in 1:3] - [kite_pos[i] => init_kite_pos[i] for i in 1:3] - [kite_vel[i] => 0 for i in 1:3] + [Q_b_w[i] => kite.orient[i] for i in 1:4] + [ω_b[i] => kite.angular_vel[i] for i in 1:3] + [kite_pos[i] => kite.pos[i] for i in 1:3] + [kite_vel[i] => kite.vel[i] for i in 1:3] ] return eqs, defaults end """ - scalar_eqs!(s, eqs, measure; R_b_w, wind_vec_gnd, va_kite_b, kite_pos, kite_vel, kite_acc) + scalar_eqs!(s, eqs, measure; R_b_w, wind_vec_gnd, va_kite_b, kite_pos, kite_vel, kite_acc, twist_angle, twist_ω) Generate equations for scalar quantities like elevation, azimuth, heading and course angles. - -# Arguments -- `s::RamAirKite`: The kite system state -- `eqs`: Current system equations -- `measure::Measurement`: Current measurement data -- `R_b_w`: Body to world rotation matrix -- `wind_vec_gnd`: Ground wind vector -- `va_kite_b`: Apparent wind velocity in body frame -- `kite_pos`: Kite position vector -- `kite_vel`: Kite velocity vector -- `kite_acc`: Kite acceleration vector - -# Returns -- Updated system equations including: - - Heading angle from x-axis - - Elevation angle - - Azimuth angle - - Course angle - - Angular velocities and accelerations -""" -function scalar_eqs!(s, eqs, measure; R_b_w, wind_vec_gnd, va_kite_b, kite_pos, kite_vel, kite_acc) + + # Arguments + - `s::RamAirKite`: The kite system state + - `eqs`: Current system equations + - `measure::Measurement`: Current measurement data + - `R_b_w`: Body to world rotation matrix + - `wind_vec_gnd`: Ground wind vector + - `va_kite_b`: Apparent wind velocity in body frame + - `kite_pos`: Kite position vector + - `kite_vel`: Kite velocity vector + - `kite_acc`: Kite acceleration vector + + # Returns + - Updated system equations including: + - Heading angle from x-axis + - Elevation angle + - Azimuth angle + - Course angle + - Angular velocities and accelerations + """ +function scalar_eqs!(s, eqs, measure; R_b_w, wind_vec_gnd, va_kite_b, kite_pos, kite_vel, kite_acc, twist_angle, twist_ω, ω_b, α_b) @parameters wind_scale_gnd = s.set.v_wind @parameters measured_wind_dir_gnd = measure.wind_dir_gnd @variables begin @@ -641,16 +654,18 @@ function scalar_eqs!(s, eqs, measure; R_b_w, wind_vec_gnd, va_kite_b, kite_pos, end eqs = [ eqs - e_x ~ R_b_w * [1, 0, 0] - e_y ~ R_b_w * [0, 1, 0] - e_z ~ R_b_w * [0, 0, 1] + e_x ~ R_b_w[:,1] + e_y ~ R_b_w[:,2] + e_z ~ R_b_w[:,3] wind_vec_gnd ~ wind_scale_gnd * rotate_around_z([1, 0, 0], measured_wind_dir_gnd) - wind_vel_kite ~ AtmosphericModels.calc_wind_factor(s.am, kite_pos[3], s.set.profile_law) * wind_vec_gnd + wind_vel_kite ~ AtmosphericModels.calc_wind_factor(s.am, kite_pos[3], s.set.profile_law) * wind_vec_gnd va_kite ~ wind_vel_kite - kite_vel va_kite_b ~ R_b_w' * va_kite ] @variables begin - heading_y(t) + heading_x(t) + turn_rate(t)[1:3] + turn_acc(t)[1:3] azimuth(t) azimuth_vel(t) azimuth_acc(t) @@ -663,15 +678,38 @@ function scalar_eqs!(s, eqs, measure; R_b_w, wind_vec_gnd, va_kite_b, kite_pos, sphere_pos(t)[1:2, 1:2] sphere_vel(t)[1:2, 1:2] sphere_acc(t)[1:2, 1:2] + angle_of_attack(t) + simple_twist_angle(t)[1:2] + simple_twist_ω(t)[1:2] + R_v_w(t)[1:3, 1:3] + view_z(t)[1:3] + view_y(t)[1:3] + view_x(t)[1:3] + distance(t) + distance_vel(t) + distance_acc(t) end x, y, z = kite_pos x´, y´, z´ = kite_vel x´´, y´´, z´´ = kite_acc + half_len = length(twist_angle)÷2 + eqs = [ eqs - heading_y ~ atan(e_x[2]/e_x[1]) + view_z ~ sym_normalize(kite_pos) + view_y ~ view_z × e_y + view_x ~ view_y × view_z + R_v_w[:,1] ~ view_x + R_v_w[:,2] ~ view_y + R_v_w[:,3] ~ view_z + heading_x ~ calc_heading(e_x) + turn_rate ~ R_v_w' * (R_b_w * ω_b) # Project angular velocity onto view frame + turn_acc ~ R_v_w' * (R_b_w * α_b) + distance ~ norm(kite_pos) + distance_vel ~ kite_vel ⋅ view_z + distance_acc ~ kite_acc ⋅ view_z elevation ~ atan(z / x) # elevation_vel = d/dt(atan(z/x)) = (x*ż' - z*ẋ')/(x^2 + z^2) according to wolframalpha @@ -686,12 +724,18 @@ function scalar_eqs!(s, eqs, measure; R_b_w, wind_vec_gnd, va_kite_b, kite_pos, x_acc ~ kite_acc ⋅ e_x y_acc ~ kite_acc ⋅ e_y course ~ atan(-azimuth_vel, elevation_vel) + + angle_of_attack ~ calc_angle_of_attack(va_kite_b) + 0.5twist_angle[half_len] + 0.5twist_angle[half_len+1] + simple_twist_angle[1] ~ sum(twist_angle[1:half_len]) / half_len + simple_twist_angle[2] ~ sum(twist_angle[half_len+1:end]) / half_len + simple_twist_ω[1] ~ sum(twist_ω[1:half_len]) / half_len + simple_twist_ω[2] ~ sum(twist_ω[half_len+1:end]) / half_len ] return eqs end """ - linear_vsm_eqs!(s, eqs; aero_force_b, aero_moment_b, group_aero_moment, init_va_b, twist_angle, va_kite_b, ω_b) +linear_vsm_eqs!(s, eqs; aero_force_b, aero_moment_b, group_aero_moment, init_va_b, twist_angle, va_kite_b, ω_b) Generate linearized aerodynamic equations using the Vortex Step Method (VSM). Uses linearization around current operating point to approximate aerodynamic forces @@ -710,23 +754,23 @@ and moments. The Jacobian is computed using the VSM solver. # Returns - Updated system equations including linearized aerodynamics: - - Force and moment calculations - - Group moment distributions - - Jacobian matrix for state derivatives +- Force and moment calculations +- Group moment distributions +- Jacobian matrix for state derivatives """ -function linear_vsm_eqs!(s, eqs; aero_force_b, aero_moment_b, group_aero_moment, init_va_b, twist_angle, va_kite_b, ω_b) +function linear_vsm_eqs!(s, eqs, guesses; aero_force_b, aero_moment_b, group_aero_moment, init_va_b, twist_angle, va_kite_b, ω_b) sol = s.vsm_solver.sol @assert length(s.point_system.groups) == length(sol.group_moment_dist) - y_ = [zeros(4); init_va_b; zeros(3)] + y_ = [init_va_b; zeros(length(s.point_system.groups)); zeros(3)] jac_, x_ = VortexStepMethod.linearize( s.vsm_solver, s.aero, y_; - theta_idxs=1:4, - va_idxs=5:7, - omega_idxs=8:10, - moment_frac=s.bridle_fracs[s.point_system.groups[1].fixed_index]) + va_idxs=1:3, + omega_idxs=4:6, + theta_idxs=7:6+length(s.point_system.groups), + moment_frac=s.point_system.groups[1].moment_frac) @parameters begin last_y[eachindex(y_)] = y_ @@ -741,180 +785,36 @@ function linear_vsm_eqs!(s, eqs; aero_force_b, aero_moment_b, group_aero_moment, eqs = [ eqs - y ~ [twist_angle; va_kite_b; ω_b] + y ~ [va_kite_b; ω_b; twist_angle] dy ~ y - last_y [aero_force_b; aero_moment_b; group_aero_moment] ~ last_x + vsm_jac * dy ] - return eqs -end - -function init_unknowns_vec!( - s::RamAirKite, - system::PointMassSystem, - vec::Vector{SimFloat}, - init_Q_b_w, - init_kite_pos; - non_observed=true -) - !s.set.quasi_static && non_observed && (length(vec) != length(s.integrator.u)) && - throw(ArgumentError("Unknowns of length $(length(s.integrator.u)) but vector provided of length $(length(vec))")) - - @unpack points, groups, segments, pulleys, winches = system - vec_idx = 1 - - if non_observed - for point in points - if point.type == DYNAMIC - for i in 1:3 - vec[vec_idx] = point.pos_w[i] - vec_idx += 1 - end - for i in 1:3 # TODO: add speed to init - vec[vec_idx] = 0.0 - vec_idx += 1 - end - end - end - - for group in groups - if group.type == DYNAMIC - vec[vec_idx] = 0 - vec_idx += 1 - vec[vec_idx] = 0 - vec_idx += 1 - end - end - - for pulley in pulleys - if pulley.type == DYNAMIC - vec[vec_idx] = segments[pulley.segments[1]].l0 - vec_idx += 1 - vec[vec_idx] = 0 - vec_idx += 1 - end - end - end - - for winch in winches - vec[vec_idx] = winch.tether_length - vec_idx += 1 - vec[vec_idx] = 0 - vec_idx += 1 - end - - for i in 1:4 - vec[vec_idx] = init_Q_b_w[i] - vec_idx += 1 - end - for i in 1:3 - vec[vec_idx] = 0 - vec_idx += 1 - end - for i in 1:3 - vec[vec_idx] = init_kite_pos[i] - vec_idx += 1 + if s.set.quasi_static + guesses = [guesses; [y[i] => y_[i] for i in eachindex(y_)]] end - for i in 1:3 - vec[vec_idx] = 0 - vec_idx += 1 - end - non_observed && (vec_idx-1 != length(vec)) && - throw(ArgumentError("Unknowns vec is of length $(length(vec)) but the last index is $(vec_idx-1)")) - nothing -end -function get_unknowns(s::RamAirKite) - vec = Num[] - vec = get_stiff_unknowns(s, vec) - vec = get_nonstiff_unknowns(s, vec) - !s.set.quasi_static && (length(vec) != length(s.integrator.u)) && - throw(ArgumentError("Integrator unknowns of length $(length(s.integrator.u)) should equal vec of length $(length(vec))")) - return vec + return eqs, guesses end -function get_stiff_unknowns(s, vec=Num[]) - @unpack points, groups, segments, pulleys, winches = s.point_system - for point in points - for i in 1:3 - point.type == DYNAMIC && push!(vec, s.sys.pos[i, point.idx]) - end - for i in 1:3 # TODO: add speed to init - point.type == DYNAMIC && push!(vec, s.sys.vel[i, point.idx]) - end - end - for group in groups - group.type == DYNAMIC && push!(vec, s.sys.free_twist_angle[group.idx]) - group.type == DYNAMIC && push!(vec, s.sys.twist_ω[group.idx]) - end - for pulley in pulleys - pulley.type == DYNAMIC && push!(vec, s.sys.pulley_l0[pulley.idx]) - pulley.type == DYNAMIC && push!(vec, s.sys.pulley_vel[pulley.idx]) - end - return vec -end - -function get_nonstiff_unknowns(s, vec=Num[]) - @unpack points, groups, segments, pulleys, winches = s.point_system - for winch in winches - push!(vec, s.sys.tether_length[winch.idx]) - push!(vec, s.sys.tether_vel[winch.idx]) - end - [push!(vec, s.sys.Q_b_w[i]) for i in 1:4] - [push!(vec, s.sys.ω_b[i]) for i in 1:3] - [push!(vec, s.sys.kite_pos[i]) for i in 1:3] - [push!(vec, s.sys.kite_vel[i]) for i in 1:3] - return vec -end - - -""" - create_sys!(s::RamAirKite, system::PointMassSystem, measure::Measurement; - init_Q_b_w, init_kite_pos, init_va_b) -> Tuple - -Create the complete ODE system for the ram air kite model by combining force equations, -aerodynamic equations, differential equations and scalar equations. - -The system is built up in stages: -1. Force equations for springs, drag and pulleys -2. Linearized aerodynamic equations using VSM -3. Differential equations for kinematics and dynamics -4. Scalar equations for angles and derived quantities - -# Arguments -- `s::RamAirKite`: The kite system state object -- `system::PointMassSystem`: The point mass representation of the kite system -- `measure::Measurement`: Initial measurement data -- `init_Q_b_w`: Initial quaternion orientation of kite body frame relative to world frame -- `init_kite_pos`: Initial kite position vector in world frame -- `init_va_b`: Initial apparent wind velocity vector in body frame - -# Returns -Tuple containing: -- `sys::ODESystem`: The complete ModelingToolkit ODE system -- `defaults::Vector{Pair{Num,Real}}`: Default values for state variables -- `guesses::Vector{Pair{Num,Real}}`: Initial guesses for solver -""" -function create_sys!(s::RamAirKite, system::PointMassSystem, measure::Measurement; init_Q_b_w, init_kite_pos, init_va_b) +function create_sys!(s::RamAirKite, system::PointMassSystem, measure::Measurement; init_va_b) eqs = [] defaults = Pair{Num, Real}[] guesses = Pair{Num, Real}[] + init_set_values = measure.set_values @parameters begin - steady = false - # measured_sphere_pos[1:2, 1:2] = measure.sphere_pos - # measured_sphere_vel[1:2, 1:2] = measure.sphere_vel - # measured_sphere_acc[1:2, 1:2] = measure.sphere_acc - # measured_tether_length[1:3] = measure.tether_length - # measured_tether_vel[1:3] = measure.tether_vel - # measured_tether_acc[1:3] = measure.tether_acc + stabilize = false + fix_nonstiff = false end @variables begin # potential differential variables + set_values(t)[eachindex(system.winches)] = init_set_values kite_pos(t)[1:3] # xyz pos of kite in world frame kite_vel(t)[1:3] kite_acc(t)[1:3] ω_b(t)[1:3] # turn rate in principal frame + α_b(t)[1:3] # rotations and frames R_b_w(t)[1:3, 1:3] # rotation of the kite body frame relative to the world frame @@ -923,6 +823,7 @@ function create_sys!(s::RamAirKite, system::PointMassSystem, measure::Measuremen aero_force_b(t)[1:3] aero_moment_b(t)[1:3] twist_angle(t)[eachindex(system.groups)] + twist_ω(t)[eachindex(system.groups)] group_aero_moment(t)[eachindex(system.groups)] wind_vec_gnd(t)[1:3] va_kite_b(t)[1:3] @@ -930,11 +831,11 @@ function create_sys!(s::RamAirKite, system::PointMassSystem, measure::Measuremen eqs, defaults, guesses, tether_kite_force, tether_kite_moment = force_eqs!(s, system, eqs, defaults, guesses; - R_b_w, kite_pos, kite_vel, wind_vec_gnd, group_aero_moment, twist_angle, steady) - eqs = linear_vsm_eqs!(s, eqs; aero_force_b, aero_moment_b, group_aero_moment, init_va_b, twist_angle, va_kite_b, ω_b) + R_b_w, kite_pos, kite_vel, wind_vec_gnd, group_aero_moment, twist_angle, twist_ω, stabilize, set_values, fix_nonstiff) + eqs, guesses = linear_vsm_eqs!(s, eqs, guesses; aero_force_b, aero_moment_b, group_aero_moment, init_va_b, twist_angle, va_kite_b, ω_b) eqs, defaults = diff_eqs!(s, eqs, defaults; tether_kite_force, tether_kite_moment, aero_force_b, aero_moment_b, - ω_b, R_b_w, kite_pos, kite_vel, kite_acc, init_Q_b_w, init_kite_pos, steady) - eqs = scalar_eqs!(s, eqs, measure; R_b_w, wind_vec_gnd, va_kite_b, kite_pos, kite_vel, kite_acc) + ω_b, α_b, R_b_w, kite_pos, kite_vel, kite_acc, stabilize, fix_nonstiff) + eqs = scalar_eqs!(s, eqs, measure; R_b_w, wind_vec_gnd, va_kite_b, kite_pos, kite_vel, kite_acc, twist_angle, twist_ω, ω_b, α_b) # te_I = (1/3 * (s.set.mass/8) * te_length^2) # # -damping / I * ω = α_damping @@ -960,6 +861,14 @@ function create_sys!(s::RamAirKite, system::PointMassSystem, measure::Measuremen @info "Creating ODESystem" # @named sys = ODESystem(eqs, t; discrete_events) @time @named sys = ODESystem(eqs, t) - return sys, defaults, guesses -end + defaults = [ + defaults + [set_values[i] => [-50.0, -1.0, -1.0][i] for i in 1:3] + ] + + s.defaults = defaults + s.guesses = guesses + s.full_sys = sys + return set_values +end diff --git a/src/point_mass_system.jl b/src/point_mass_system.jl index 0c5d35f4c..74025e027 100644 --- a/src/point_mass_system.jl +++ b/src/point_mass_system.jl @@ -1,18 +1,287 @@ -# Copyright (c) 2024 Bart van de Lint -# SPDX-License-Identifier: MIT +# Copyright (c) 2025 Bart van de Lint +# SPDX-License-Identifier: MPL-2.0 -function VortexStepMethod.RamAirWing(set::Settings; prn=true) +function VortexStepMethod.RamAirWing(set::Settings; prn=true, kwargs...) obj_path = joinpath(dirname(get_data_path()), set.model) dat_path = joinpath(dirname(get_data_path()), set.foil_file) - return RamAirWing(obj_path, dat_path; mass=set.mass, crease_frac=set.crease_frac, align_to_principal=true, prn) + return RamAirWing(obj_path, dat_path; + mass=set.mass, crease_frac=set.crease_frac, n_groups=4, + align_to_principal=true, prn, kwargs... + ) end +""" + SegmentType `POWER` `STEERING` `BRIDLE` + +Type of segment. + +# Elements +- POWER: Belongs to a power line +- STEERING: Belongs to a steering line +- BRIDLE: Belongs to the bridle +""" +@enum SegmentType begin + POWER + STEERING + BRIDLE +end + +""" + DynamicsType `DYNAMIC` `STATIC` `KITE` `WINCH` + +Enumeration of the models that are attached to a point. + +# Elements +- DYNAMIC: Belongs to a dynamic tether model +- STATIC: Belongs to a static tether model +- KITE: Rigid body +- WINCH: Winch +""" +@enum DynamicsType begin + DYNAMIC + STATIC + KITE + WINCH +end + +""" + mutable struct Point + +A normal freely moving tether point. + +$(TYPEDFIELDS) +""" +mutable struct Point + idx::Int16 + pos_b::KVec3 # pos relative to kite COM in body frame + pos_w::KVec3 # pos in world frame + vel_w::KVec3 # vel in world frame + type::DynamicsType +end +function Point(idx, pos_b, type, vel_w=zeros(KVec3)) + Point(idx, pos_b, copy(pos_b), vel_w, type) +end + +""" + struct KitePointGroup + +Set of bridle lines that share the same twist angle and trailing edge angle. + +$(TYPEDFIELDS) +""" +mutable struct KitePointGroup + idx::Int16 + points::Vector{Int16} + le_pos::KVec3 # point which the group rotates around under kite deformation + chord::KVec3 # chord vector in body frame which the group rotates around under kite deformation + y_airf::KVec3 # spanwise vector in local panel frame which the group rotates around under kite deformation + type::DynamicsType + moment_frac::SimFloat + twist::SimFloat + twist_vel::SimFloat +end +function KitePointGroup(idx, points, wing, gamma, type, moment_frac) + le_pos = [wing.le_interp[i](gamma) for i in 1:3] + chord = [wing.te_interp[i](gamma) for i in 1:3] .- le_pos + y_airf = normalize([wing.le_interp[i](gamma-0.01) for i in 1:3] - le_pos) + KitePointGroup(idx, points, le_pos, chord, y_airf, type, moment_frac, 0.0, 0.0) +end + +""" + mutable struct Segment + +A segment from one point index to another point index. + +$(TYPEDFIELDS) +""" +mutable struct Segment + idx::Int16 + points::Tuple{Int16, Int16} + type::SegmentType + l0::SimFloat + compression_frac::SimFloat + diameter::SimFloat +end +function Segment(idx, points, type, l0=zero(SimFloat), compression_frac=0.1) + Segment(idx, points, type, l0, compression_frac, zero(SimFloat)) +end + +""" + mutable struct Pulley + +A pulley described by two segments with the common point of the segments being the pulley. + +$(TYPEDFIELDS) +""" +mutable struct Pulley + idx::Int16 + segments::Tuple{Int16, Int16} + type::DynamicsType + sum_length::SimFloat + length::SimFloat + vel::SimFloat + function Pulley(idx, segments, type) + new(idx, segments, type, 0.0, 0.0, 0.0) + end +end + +""" + struct Tether + +A set of segments making a flexible tether. The winch point should only be part of one segment. + +$(TYPEDFIELDS) +""" +struct Tether + idx::Int16 + segments::Vector{Int16} + winch_point::Int16 +end + +""" + mutable struct Winch + +A set of tethers or just one tether connected to a winch. + +$(TYPEDFIELDS) +""" +mutable struct Winch + idx::Int16 + model::AbstractWinchModel + tethers::Vector{Int16} + tether_length::SimFloat + tether_vel::SimFloat + function Winch(idx, model, tethers, tether_length, tether_vel=0.0) + new(idx, model, tethers, tether_length, tether_vel) + end +end + +@with_kw struct Kite + orient::KVec4 = zeros(KVec4) + angular_vel::KVec3 = zeros(KVec3) + pos::KVec3 = zeros(KVec3) + vel::KVec3 = zeros(KVec3) +end + +""" + struct PointMassSystem + +A discrete mass-spring-damper representation of a kite system, where point masses +connected by elastic segments model the kite and tether dynamics: + +- `points::Vector{Point}`: Point masses representing: + - Kite attachment points + - Dynamic bridle/tether points + - Fixed ground anchor points +- `groups::Vector{KitePointGroup}`: Collections of points that move together, + according to kite deformation (twist and trailing edge deflection) +- `segments::Vector{Segment}`: Spring-damper elements between points +- `pulleys::Vector{Pulley}`: Elements that redistribute line lengths +- `tethers::Vector{Tether}`: Groups of segments with a common unstretched length +- `winches::Vector{Winch}`: Ground-based winches that control the tether lengths + +See also: [`Point`](@ref), [`Segment`](@ref), [`KitePointGroup`](@ref), [`Pulley`](@ref) +""" +struct PointMassSystem + name::String + points::Vector{Point} + groups::Vector{KitePointGroup} + segments::Vector{Segment} + pulleys::Vector{Pulley} + tethers::Vector{Tether} + winches::Vector{Winch} + kite::Kite + function PointMassSystem(name, points, groups, segments, pulleys, tethers, winches, kite=Kite()) + for (i, point) in enumerate(points) + @assert point.idx == i + end + for (i, group) in enumerate(groups) + @assert group.idx == i + end + for (i, segment) in enumerate(segments) + @assert segment.idx == i + end + for (i, pulley) in enumerate(pulleys) + @assert pulley.idx == i + end + for (i, tether) in enumerate(tethers) + @assert tether.idx == i + end + for (i, winch) in enumerate(winches) + @assert winch.idx == i + end + return new(name, points, groups, segments, pulleys, tethers, winches, kite) + end +end + + function PointMassSystem(set::Settings, wing::RamAirWing) - # TODO: move as much of the code as possible from create_point_mass_system to other places, to make model creation easier. - # 1. move bridle gamma calculation - # 2. ... + length(set.bridle_fracs) != 4 && throw(ArgumentError("4 bridle fracs should be provided for all models.")) + + if set.physical_model == "ram" + return create_ram_point_system(set, wing) + elseif set.physical_model == "simple_ram" + return create_simple_ram_point_system(set, wing) + else + throw(ArgumentError("Undefined physical model")) + end +end + +function calc_pos(wing::RamAirWing, gamma, frac) + le_pos = [wing.le_interp[i](gamma) for i in 1:3] + chord = [wing.te_interp[i](gamma) for i in 1:3] .- le_pos + pos = le_pos .+ chord .* frac + return pos +end +function create_tether(idx, set, points, segments, tethers, attach_point, type, dynamics_type) + segment_idxs = Int16[] + winch_pos = find_z_axis_point(attach_point.pos_b, set.l_tether) + dir = winch_pos - attach_point.pos_b + for i in 1:set.segments + frac = i / set.segments + pos = attach_point.pos_b + frac * dir + i_pnt = length(points) # last point idx + i_seg = length(segments) # last segment idx + if i == 1 + points = [points; Point(1+i_pnt, pos, dynamics_type)] + segments = [segments; Segment(1+i_seg, (attach_point.idx, 1+i_pnt), type)] + elseif i == set.segments + points = [points; Point(1+i_pnt, pos, WINCH)] + segments = [segments; Segment(1+i_seg, (i_pnt, 1+i_pnt), type)] + else + points = [points; Point(1+i_pnt, pos, dynamics_type)] + segments = [segments; Segment(1+i_seg, (i_pnt, 1+i_pnt), type)] + end + push!(segment_idxs, 1+i_seg) + i_pnt = length(points) + end + winch_point_idx = points[end].idx + tethers = [tethers; Tether(idx, segment_idxs, winch_point_idx)] + return points, segments, tethers, tethers[end].idx +end + +function cad_to_body_frame(wing::RamAirWing, pos) + return wing.R_cad_body * (pos + wing.T_cad_body) +end + +""" +Find the point on the z-axis with distance l from P in the negative direction +""" +function find_z_axis_point(P, l) + # Distance formula: sqrt(xₐ² + yₐ² + (z - zₐ)²) = l + # Square both sides: xₐ² + yₐ² + (z - zₐ)² = l² + # Solve for z: (z - zₐ)² = l² - xₐ² - yₐ² + # z - zₐ = ±sqrt(l² - xₐ² - yₐ²) + # z = zₐ ± sqrt(l² - xₐ² - yₐ²) + d = l^2 - P[1]^2 - P[2]^2 + d < 0 && error("No real solution: l is too small") + z = P[3] - sqrt(d) + return [0, 0, z] +end + +function create_ram_point_system(set::Settings, wing::RamAirWing) points = Point[] groups = KitePointGroup[] segments = Segment[] @@ -22,7 +291,7 @@ function PointMassSystem(set::Settings, wing::RamAirWing) attach_points = Point[] - bridle_top_left = [wing.R_cad_body * (set.top_bridle_points[i] + wing.T_cad_body) for i in eachindex(set.top_bridle_points)] # cad to kite frame + bridle_top_left = [cad_to_body_frame(wing, set.top_bridle_points[i]) for i in eachindex(set.top_bridle_points)] bridle_top_right = [bridle_top_left[i] .* [1, -1, 1] for i in eachindex(set.top_bridle_points)] dynamics_type = set.quasi_static ? STATIC : DYNAMIC @@ -31,124 +300,186 @@ function PointMassSystem(set::Settings, wing::RamAirWing) i_pnt = length(points) # last point idx i_seg = length(segments) # last segment idx i_pul = length(pulleys) # last pulley idx + i_grp = length(groups) # last group idx - i = 1 - for gamma in gammas # 2 gammas per bridle system - le_pos = [wing.le_interp[i](gamma) for i in 1:3] - chord = [wing.te_interp[i](gamma) for i in 1:3] .- le_pos - y_airf = normalize([wing.le_interp[i](gamma-0.01) for i in 1:3] - le_pos) - - point_idxs = Int16[] - for frac in set.bridle_fracs # 4 fracs - pos = le_pos .+ chord .* frac - points = [points; Point(i+i_pnt, pos, KITE)] - push!(point_idxs, points[end].idx) - i += 1 - end - - i_grp = 1 + length(groups) - groups = [groups; KitePointGroup(i_grp, point_idxs, 1, chord, y_airf, DYNAMIC)] - end + # ==================== CREATE DEFORMING KITE GROUPS ==================== # + points = [ + points + Point(1+i_pnt, calc_pos(wing, gammas[1], set.bridle_fracs[1]), KITE) + Point(2+i_pnt, calc_pos(wing, gammas[1], set.bridle_fracs[3]), KITE) + Point(3+i_pnt, calc_pos(wing, gammas[1], set.bridle_fracs[4]), KITE) + + Point(4+i_pnt, calc_pos(wing, gammas[2], set.bridle_fracs[1]), KITE) + Point(5+i_pnt, calc_pos(wing, gammas[2], set.bridle_fracs[3]), KITE) + Point(6+i_pnt, calc_pos(wing, gammas[2], set.bridle_fracs[4]), KITE) + ] + groups = [ + groups + KitePointGroup(1+i_grp, [1+i_pnt, 2+i_pnt, 3+i_pnt], wing, gammas[1], DYNAMIC, set.bridle_fracs[2]) + KitePointGroup(2+i_grp, [4+i_pnt, 5+i_pnt, 6+i_pnt], wing, gammas[2], DYNAMIC, set.bridle_fracs[2]) + ] + # ==================== CREATE PULLEY BRIDLE SYSTEM ==================== # + # TODO: add initial rotation around y-axis points = [ points - Point(9+i_pnt, bridle_top[1], dynamics_type) - Point(10+i_pnt, bridle_top[2], dynamics_type) - Point(11+i_pnt, bridle_top[3], dynamics_type) - Point(12+i_pnt, bridle_top[4], dynamics_type) + Point(7+i_pnt, bridle_top[1], dynamics_type) + Point(8+i_pnt, bridle_top[2], KITE) + Point(9+i_pnt, bridle_top[3], dynamics_type) + Point(10+i_pnt, bridle_top[4], dynamics_type) - Point(13+i_pnt, bridle_top[2] .+ [0, 0, -1], dynamics_type) + Point(11+i_pnt, bridle_top[2] .+ [0, 0, -1], dynamics_type) - Point(14+i_pnt, bridle_top[1] .+ [0, 0, -2], dynamics_type) - Point(15+i_pnt, bridle_top[3] .+ [0, 0, -2], dynamics_type) + Point(12+i_pnt, bridle_top[1] .+ [0, 0, -2], dynamics_type) + Point(13+i_pnt, bridle_top[3] .+ [0, 0, -2], dynamics_type) - Point(16+i_pnt, bridle_top[1] .+ [0, 0, -3], dynamics_type) - Point(17+i_pnt, bridle_top[3] .+ [0, 0, -3], dynamics_type) + Point(14+i_pnt, bridle_top[1] .+ [0, 0, -4], dynamics_type) + Point(15+i_pnt, bridle_top[3] .+ [0, 0, -4], dynamics_type) ] - l1 = norm(points[9+i_pnt].pos_b - points[1+i_pnt].pos_b) - l2 = norm(points[9+i_pnt].pos_b - points[5+i_pnt].pos_b) segments = [ segments - Segment(1+i_seg, (1+i_pnt, 9+i_pnt), BRIDLE, l1) - Segment(2+i_seg, (2+i_pnt, 10+i_pnt), BRIDLE, l1) - Segment(3+i_seg, (3+i_pnt, 11+i_pnt), BRIDLE, l1) - Segment(4+i_seg, (4+i_pnt, 12+i_pnt), BRIDLE, l1) - - Segment(5+i_seg, (5+i_pnt, 9+i_pnt), BRIDLE, l2) - Segment(6+i_seg, (6+i_pnt, 10+i_pnt), BRIDLE, l2) - Segment(7+i_seg, (7+i_pnt, 11+i_pnt), BRIDLE, l2) - Segment(8+i_seg, (8+i_pnt, 12+i_pnt), BRIDLE, l2) - - Segment(9+i_seg, (9+i_pnt, 14+i_pnt), BRIDLE, 2) - Segment(10+i_seg, (10+i_pnt, 13+i_pnt), BRIDLE, 1) - Segment(11+i_seg, (11+i_pnt, 15+i_pnt), BRIDLE, 2) - Segment(12+i_seg, (12+i_pnt, 17+i_pnt), BRIDLE, 3) + Segment(1+i_seg, (1+i_pnt, 7+i_pnt), BRIDLE) + Segment(2+i_seg, (2+i_pnt, 9+i_pnt), BRIDLE) + Segment(3+i_seg, (3+i_pnt, 10+i_pnt), BRIDLE) + + Segment(4+i_seg, (4+i_pnt, 7+i_pnt), BRIDLE) + Segment(5+i_seg, (5+i_pnt, 9+i_pnt), BRIDLE) + Segment(6+i_seg, (6+i_pnt, 10+i_pnt), BRIDLE) + + Segment(7+i_seg, (7+i_pnt, 12+i_pnt), BRIDLE, 2) + Segment(8+i_seg, (8+i_pnt, 11+i_pnt), BRIDLE, 1) + Segment(9+i_seg, (9+i_pnt, 13+i_pnt), BRIDLE, 2) + Segment(10+i_seg, (10+i_pnt, 15+i_pnt), BRIDLE, 4) - Segment(13+i_seg, (13+i_pnt, 14+i_pnt), BRIDLE, 1) - Segment(14+i_seg, (13+i_pnt, 15+i_pnt), BRIDLE, 1) + Segment(11+i_seg, (11+i_pnt, 12+i_pnt), BRIDLE, 1) + Segment(12+i_seg, (11+i_pnt, 13+i_pnt), BRIDLE, 1) - Segment(15+i_seg, (14+i_pnt, 16+i_pnt), BRIDLE, 1) - Segment(16+i_seg, (15+i_pnt, 16+i_pnt), BRIDLE, 1) - Segment(17+i_seg, (15+i_pnt, 17+i_pnt), BRIDLE, 1) + Segment(13+i_seg, (12+i_pnt, 14+i_pnt), BRIDLE, 2) + Segment(14+i_seg, (13+i_pnt, 14+i_pnt), BRIDLE, 2) + Segment(15+i_seg, (13+i_pnt, 15+i_pnt), BRIDLE, 2) ] pulleys = [ pulleys - Pulley(1+i_pul, (13+i_seg, 14+i_seg), dynamics_type) - Pulley(2+i_pul, (16+i_seg, 17+i_seg), dynamics_type) + Pulley(1+i_pul, (11+i_seg, 12+i_seg), dynamics_type) + Pulley(2+i_pul, (14+i_seg, 15+i_seg), dynamics_type) ] push!(attach_points, points[end-1]) push!(attach_points, points[end]) return nothing end - function create_tether(attach_point, type) - l0 = set.l_tether / set.segments - segment_idxs = Int16[] - for i in 1:set.segments - frac = i / set.segments - pos = [(1-frac) * attach_point.pos_b[1], - (1-frac) * attach_point.pos_b[2], - attach_point.pos_b[3] - i*l0] - i_pnt = length(points) # last point idx - i_seg = length(segments) # last segment idx - if i == 1 - points = [points; Point(1+i_pnt, pos, dynamics_type)] - segments = [segments; Segment(1+i_seg, (attach_point.idx, 1+i_pnt), type)] - elseif i == set.segments - points = [points; Point(1+i_pnt, pos, WINCH)] - segments = [segments; Segment(1+i_seg, (i_pnt, 1+i_pnt), type)] - else - points = [points; Point(1+i_pnt, pos, dynamics_type)] - segments = [segments; Segment(1+i_seg, (i_pnt, 1+i_pnt), type)] - end - push!(segment_idxs, 1+i_seg) - i_pnt = length(points) - end - i_tether = length(tethers) - winch_point_idx = points[end].idx - tethers = [tethers; Tether(1+i_tether, segment_idxs, winch_point_idx)] - return tethers[end].idx - end + # function create_tether(attach_point, type) + # l0 = set.l_tether / set.segments + # segment_idxs = Int16[] + # winch_pos = find_z_axis_point(attach_point.pos_b, set.l_tether) + # dir = winch_pos - attach_point.pos_b + # for i in 1:set.segments + # frac = i / set.segments + # pos = attach_point.pos_b + frac * dir + # i_pnt = length(points) # last point idx + # i_seg = length(segments) # last segment idx + # if i == 1 + # points = [points; Point(1+i_pnt, pos, dynamics_type)] + # segments = [segments; Segment(1+i_seg, (attach_point.idx, 1+i_pnt), type)] + # elseif i == set.segments + # points = [points; Point(1+i_pnt, pos, WINCH)] + # segments = [segments; Segment(1+i_seg, (i_pnt, 1+i_pnt), type)] + # else + # points = [points; Point(1+i_pnt, pos, dynamics_type)] + # segments = [segments; Segment(1+i_seg, (i_pnt, 1+i_pnt), type)] + # end + # push!(segment_idxs, 1+i_seg) + # i_pnt = length(points) + # end + # i_tether = length(tethers) + # winch_point_idx = points[end].idx + # tethers = [tethers; Tether(1+i_tether, segment_idxs, winch_point_idx)] + # return tethers[end].idx + # end + + gammas = [-3/4, -1/4, 1/4, 3/4] * wing.gamma_tip + create_bridle(bridle_top_left, gammas[[1,2]]) + create_bridle(bridle_top_right, gammas[[4,3]]) + + + # left_power_idx = create_tether(attach_points[1], POWER) + # right_power_idx = create_tether(attach_points[3], POWER) + # left_steering_idx = create_tether(attach_points[2], STEERING) + # right_steering_idx = create_tether(attach_points[4], STEERING) + + points, segments, tethers, left_power_idx = create_tether(1, set, points, segments, tethers, attach_points[1], POWER, dynamics_type) + points, segments, tethers, right_power_idx = create_tether(2, set, points, segments, tethers, attach_points[3], POWER, dynamics_type) + points, segments, tethers, left_steering_idx = create_tether(3, set, points, segments, tethers, attach_points[2], STEERING, dynamics_type) + points, segments, tethers, right_steering_idx = create_tether(4, set, points, segments, tethers, attach_points[4], STEERING, dynamics_type) + + winches = [winches; Winch(1, TorqueControlledMachine(set), [left_power_idx, right_power_idx], set.l_tether)] + winches = [winches; Winch(2, TorqueControlledMachine(set), [left_steering_idx], set.l_tether)] + winches = [winches; Winch(3, TorqueControlledMachine(set), [right_steering_idx], set.l_tether)] + return PointMassSystem(set.physical_model, points, groups, segments, pulleys, tethers, winches) +end + +function create_simple_ram_point_system(set::Settings, wing::RamAirWing) + points = Point[] + groups = KitePointGroup[] + segments = Segment[] + pulleys = Pulley[] + tethers = Tether[] + winches = Winch[] + + dynamics_type = set.quasi_static ? STATIC : DYNAMIC gammas = [-3/4, -1/4, 1/4, 3/4] * wing.gamma_tip - create_bridle(bridle_top_left, gammas[1:2]) - create_bridle(bridle_top_right, gammas[3:4]) + + bridle_top_left = [wing.R_cad_body * (set.top_bridle_points[i] + wing.T_cad_body) for i in eachindex(set.top_bridle_points)] # cad to kite frame + bridle_top_right = [bridle_top_left[i] .* [1, -1, 1] for i in eachindex(set.top_bridle_points)] + + # ==================== CREATE DEFORMING KITE GROUPS ==================== # + points = [ + points + Point(1, calc_pos(wing, gammas[1], set.bridle_fracs[4]), KITE) + Point(2, calc_pos(wing, gammas[2], set.bridle_fracs[4]), KITE) + Point(3, calc_pos(wing, gammas[3], set.bridle_fracs[4]), KITE) + Point(4, calc_pos(wing, gammas[4], set.bridle_fracs[4]), KITE) + ] + groups = [ + groups + KitePointGroup(1, [1], wing, gammas[1], DYNAMIC, set.bridle_fracs[2]) + KitePointGroup(2, [2], wing, gammas[2], DYNAMIC, set.bridle_fracs[2]) + KitePointGroup(3, [3], wing, gammas[3], DYNAMIC, set.bridle_fracs[2]) + KitePointGroup(4, [4], wing, gammas[4], DYNAMIC, set.bridle_fracs[2]) + ] + # ==================== CREATE PULLEY BRIDLE SYSTEM ==================== # + points = [ + points + Point(5, bridle_top_left[2], KITE) + Point(6, bridle_top_left[4], dynamics_type) + Point(7, bridle_top_right[2], KITE) + Point(8, bridle_top_right[4], dynamics_type) + ] - left_power_idx = create_tether(attach_points[1], POWER) - right_power_idx = create_tether(attach_points[3], POWER) - left_steering_idx = create_tether(attach_points[2], STEERING) - right_steering_idx = create_tether(attach_points[4], STEERING) + segments = [ + segments + Segment(1, (1, 6), BRIDLE) + Segment(2, (2, 6), BRIDLE) + Segment(3, (3, 8), BRIDLE) + Segment(4, (4, 8), BRIDLE) + ] - winches = [winches; Winch(1, TorqueControlledMachine(set), [left_power_idx, right_power_idx])] - winches = [winches; Winch(2, TorqueControlledMachine(set), [left_steering_idx])] - winches = [winches; Winch(3, TorqueControlledMachine(set), [right_steering_idx])] + points, segments, tethers, left_power_idx = create_tether(1, set, points, segments, tethers, points[5], POWER, dynamics_type) + points, segments, tethers, right_power_idx = create_tether(2, set, points, segments, tethers, points[7], POWER, dynamics_type) + points, segments, tethers, left_steering_idx = create_tether(3, set, points, segments, tethers, points[6], STEERING, dynamics_type) + points, segments, tethers, right_steering_idx = create_tether(4, set, points, segments, tethers, points[8], STEERING, dynamics_type) - return PointMassSystem(points, groups, segments, pulleys, tethers, winches) + winches = [winches; Winch(1, TorqueControlledMachine(set), [left_power_idx, right_power_idx], set.l_tether)] + winches = [winches; Winch(2, TorqueControlledMachine(set), [left_steering_idx], set.l_tether)] + winches = [winches; Winch(3, TorqueControlledMachine(set), [right_steering_idx], set.l_tether)] + + return PointMassSystem(set.physical_model, points, groups, segments, pulleys, tethers, winches) end -function init!(system::PointMassSystem, set::Settings, R_b_w) - @unpack points, groups, segments, pulleys, tethers, winches = system +function init!(system::PointMassSystem, set::Settings, R_b_w, Q_b_w) + @unpack points, groups, segments, pulleys, tethers, winches, kite = system for segment in segments (segment.type === BRIDLE) && (segment.diameter = 0.001set.bridle_tether_diameter) @@ -162,20 +493,24 @@ function init!(system::PointMassSystem, set::Settings, R_b_w) for pulley in pulleys segment1, segment2 = segments[pulley.segments[1]], segments[pulley.segments[2]] pulley.sum_length = segment1.l0 + segment2.l0 + pulley.length = segment1.l0 + pulley.vel = 0.0 @assert !(pulley.sum_length ≈ 0) end for winch in winches - tether_length = 0.0 - for tether in tethers[winch.tethers] - for segment in segments[tether.segments] - tether_length += segment.l0 / length(winch.tethers) * 0.9999 - end - end - winch.tether_length = tether_length + winch.tether_length = set.l_tether + winch.tether_vel = 0.0 @assert !(winch.tether_length ≈ 0) end + first_moment_frac = groups[1].moment_frac + for group in groups + group.twist = 0.0 + group.twist_vel = 0.0 + @assert group.moment_frac ≈ first_moment_frac "All group.moment_frac must be the same." + end + min_z = Inf for point in points if point.pos_b[3] < min_z @@ -184,12 +519,59 @@ function init!(system::PointMassSystem, set::Settings, R_b_w) end for point in points point.pos_w .= R_b_w * [point.pos_b[1], point.pos_b[2], point.pos_b[3] - min_z] + point.vel_w .= 0.0 + end + kite.pos .= R_b_w * [0.0, 0.0, -min_z] + kite.orient .= Q_b_w + kite.vel .= 0.0 + kite.angular_vel .= 0.0 + return nothing +end + +# function update!(system::PointMassSystem, measure::Measurement) +# @unpack points, groups, segments, pulleys, tethers, winches, kite = system + +# kite.pos .= R_b_w * [0.0, 0.0, -min_z] +# kite.orient .= Q_b_w +# kite.vel .= 0.0 +# kite.angular_vel .= 0.0 + +# for (winch, tether_length, tether_vel) in zip(winches, measure.tether_length, measure.tether_vel) +# winch.tether_length = tether_length +# winch.tether_vel = tether_vel +# end + +# return nothing +# end + +const MeasureFloat = Float32 + +@with_kw mutable struct Measurement + set_values::MVector{3, MeasureFloat} = [-50., -1., -1.] + tether_length::MVector{3, MeasureFloat} = [51., 51., 51.] + tether_vel::MVector{3, MeasureFloat} = zeros(MeasureFloat, 3) + tether_acc::MVector{3, MeasureFloat} = zeros(MeasureFloat, 3) + "elevation and azimuth in spherical coordinate system with columns (left, right) and rows (elevation, azimuth)" + sphere_pos::Matrix{MeasureFloat} = deg2rad.([80.0 80.0; 1.0 -1.0]) + sphere_vel::Matrix{MeasureFloat} = zeros(MeasureFloat, 2, 2) + sphere_acc::Matrix{MeasureFloat} = zeros(MeasureFloat, 2, 2) + "positive azimuth wind direction in right-handed ENU frame relative to east / x-axis" + wind_dir_gnd::MeasureFloat = zero(MeasureFloat) +end + +function Base.getproperty(m::Measurement, val::Symbol) + if val === :elevation + sphere_pos = getfield(m, :sphere_pos) + return 0.5(sphere_pos[1, 1] + sphere_pos[1, 2]) + elseif val === :azimuth + sphere_pos = getfield(m, :sphere_pos) + return 0.5(sphere_pos[2, 1] + sphere_pos[2, 2]) + else + return getfield(m, val) end - init_kite_pos = R_b_w * [0.0, 0.0, -min_z] - return init_kite_pos end -function measure_to_q(measure::Measurement) +function measure_to_q(measure::Measurement, R_cad_body=I(3)) x = [0, 0, -1] # laying flat along x axis z = [1, 0, 0] # laying flat along x axis x = rotate_around_y(x, -measure.elevation) diff --git a/src/precompile.jl b/src/precompile.jl index 9167d6044..245d8e7cb 100644 --- a/src/precompile.jl +++ b/src/precompile.jl @@ -45,12 +45,18 @@ end kps3_::KPS3 = KPS3(KCU(se("system.yaml"))) ver = "$(VERSION.major).$(VERSION.minor)_" - input_file = joinpath(path, "..", "data", "prob_dynamic_" * ver * "3_seg.bin.default.xz") - prob_file = joinpath(path, "..", "data", "prob_dynamic_" * ver * "3_seg.bin") - output_file = joinpath(prob_file * ".default") + prob_name = "prob_" * ver * "ram_dynamic_3_seg.bin" + prob_file = joinpath(path, "..", "data", prob_name) + output_file = joinpath(path, "..", "data", prob_name * ".default") + input_file = joinpath(path, "..", "data", prob_name * ".default.xz") if isfile(input_file) && ! isfile(output_file) using CodecXz decompress_binary(input_file, output_file) + @info "Decompressed $input_file to $output_file" + elseif isfile(output_file) + @info "Output file $output_file already exists, skipping decompression." + else + @error "Input file $input_file does not exist, skipping decompression." end @assert ! isnothing(kps4_.wm) @@ -66,37 +72,27 @@ end m1 = "Manifest-v1.10.toml" m2 = "Manifest-v1.10.toml.default" end + if filecmp(m1, m2) + @info "Manifest files match, using the default xz files will work!" + else + @warn "Manifest files differ, no precompilation will be done." + end + # Check if the output file exists and is the same as the input file if isfile(output_file) && filecmp(m1, m2) local set - # Simulation parameters - dt = 0.05 - total_time = 10 # Longer simulation to see oscillations - vsm_interval = 2 - steps = Int(round(total_time / dt)) - - # Steering parameters - steering_freq = 1/2 # Hz - full left-right cycle frequency - steering_magnitude = 5.0 # Magnitude of steering input [Nm] - # Initialize model set = se("system_ram.yaml") set.segments = 3 set_values = [-50, 0.0, 0.0] # Set values of the torques of the three winches. [Nm] set.quasi_static = false - - wing = RamAirWing(set; prn=false) - aero = BodyAerodynamics([wing]) - vsm_solver = Solver(aero; solver_type=NONLIN, atol=2e-8, rtol=2e-8) - point_system = PointMassSystem(set, wing) - s = RamAirKite(set, aero, vsm_solver, point_system) - + set.physical_model = "ram" + s = RamAirKite(set) measure = Measurement() - s.set.abs_tol = 1e-5 - s.set.rel_tol = 1e-4 # Initialize at elevation measure.sphere_pos .= deg2rad.([60.0 60.0; 1.0 -1.0]) KiteModels.init_sim!(s, measure; prn=false, precompile=true) + @info "Copying $output_file to $prob_file !" cp(output_file, prob_file; force=true) end nothing diff --git a/src/ram_air_kite.jl b/src/ram_air_kite.jl index e15e084e5..058ed48cc 100644 --- a/src/ram_air_kite.jl +++ b/src/ram_air_kite.jl @@ -1,218 +1,6 @@ # Copyright (c) 2024 Bart van de Lint and Uwe Fechner # SPDX-License-Identifier: MIT -const MeasureFloat = Float32 - -@with_kw mutable struct Measurement - set_values::MVector{3, MeasureFloat} = [-50., -1., -1.] - tether_length::MVector{3, MeasureFloat} = [51., 51., 51.] - tether_vel::MVector{3, MeasureFloat} = zeros(MeasureFloat, 3) - tether_acc::MVector{3, MeasureFloat} = zeros(MeasureFloat, 3) - tether_force::MVector{3, MeasureFloat} = [540., 3., 3.] - "elevation and azimuth in spherical coordinate system with columns (left, right) and rows (elevation, azimuth)" - sphere_pos::Matrix{MeasureFloat} = deg2rad.([89.0 89.0; 1.0 -1.0]) - sphere_vel::Matrix{MeasureFloat} = zeros(MeasureFloat, 2, 2) - sphere_acc::Matrix{MeasureFloat} = zeros(MeasureFloat, 2, 2) - "positive azimuth wind direction in right-handed ENU frame relative to east / x-axis" - wind_dir_gnd::MeasureFloat = zero(MeasureFloat) -end - -function Base.getproperty(m::Measurement, val::Symbol) - if val === :elevation - sphere_pos = getfield(m, :sphere_pos) - return 0.5(sphere_pos[1, 1] + sphere_pos[1, 2]) - elseif val === :azimuth - sphere_pos = getfield(m, :sphere_pos) - return 0.5(sphere_pos[2, 1] + sphere_pos[2, 2]) - else - return getfield(m, val) - end -end - -""" - SegmentType `POWER` `STEERING` `BRIDLE` - -Type of segment. - -# Elements -- POWER: Belongs to a power line -- STEERING: Belongs to a steering line -- BRIDLE: Belongs to the bridle -""" -@enum SegmentType begin - POWER - STEERING - BRIDLE -end - -""" - DynamicsType `DYNAMIC` `STATIC` `KITE` `WINCH` - -Enumeration of the models that are attached to a point. - -# Elements -- DYNAMIC: Belongs to a dynamic tether model -- STATIC: Belongs to a static tether model -- KITE: Rigid body -- WINCH: Winch -""" -@enum DynamicsType begin - DYNAMIC - STATIC - KITE - WINCH -end - -""" - mutable struct Point - -A normal freely moving tether point. - -$(TYPEDFIELDS) -""" -mutable struct Point - idx::Int16 - pos_b::KVec3 # pos relative to kite COM in body frame - pos_w::KVec3 # pos in world frame - type::DynamicsType -end -function Point(idx, pos_b, type) - Point(idx, pos_b, copy(pos_b), type) -end - -""" - struct KitePointGroup - -Set of bridle lines that share the same twist angle and trailing edge angle. - -$(TYPEDFIELDS) -""" -struct KitePointGroup - idx::Int16 - points::Vector{Int16} - fixed_index::Int16 # point which the group rotates around under kite deformation - chord::KVec3 # chord vector in body frame which the group rotates around under kite deformation - y_airf::KVec3 # spanwise vector in local panel frame which the group rotates around under kite deformation - type::DynamicsType -end - -""" - mutable struct Segment - -A segment from one point index to another point index. - -$(TYPEDFIELDS) -""" -mutable struct Segment - idx::Int16 - points::Tuple{Int16, Int16} - type::SegmentType - l0::SimFloat - diameter::SimFloat -end -function Segment(idx, points, type) - Segment(idx, points, type, zero(SimFloat), zero(SimFloat)) -end -function Segment(idx, points, type, l0) - Segment(idx, points, type, l0, zero(SimFloat)) -end - -""" - mutable struct Pulley - -A pulley described by two segments with the common point of the segments being the pulley. - -$(TYPEDFIELDS) -""" -mutable struct Pulley - idx::Int16 - segments::Tuple{Int16, Int16} - type::DynamicsType - sum_length::SimFloat - function Pulley(idx, segments, type) - new(idx, segments, type, zero(SimFloat)) - end -end - -""" - struct Tether - -A set of segments making a flexible tether. The winch point should only be part of one segment. - -$(TYPEDFIELDS) -""" -struct Tether - idx::Int16 - segments::Vector{Int16} - winch_point::Int16 -end - -""" - mutable struct Winch - -A set of tethers or just one tether connected to a winch. - -$(TYPEDFIELDS) -""" -mutable struct Winch - idx::Int16 - model::AbstractWinchModel - tethers::Vector{Int16} - tether_length::Float64 - function Winch(idx, model, tethers) - new(idx, model, tethers, zero(Float64)) - end -end - -""" - struct PointMassSystem - -A discrete mass-spring-damper representation of a kite system, where point masses -connected by elastic segments model the kite and tether dynamics: - -- `points::Vector{Point}`: Point masses representing: - - Kite attachment points - - Dynamic bridle/tether points - - Fixed ground anchor points -- `groups::Vector{KitePointGroup}`: Collections of points that move together, - according to kite deformation (twist and trailing edge deflection) -- `segments::Vector{Segment}`: Spring-damper elements between points -- `pulleys::Vector{Pulley}`: Elements that redistribute line lengths -- `tethers::Vector{Tether}`: Groups of segments with a common unstretched length -- `winches::Vector{Winch}`: Ground-based winches that control the tether lengths - -See also: [`Point`](@ref), [`Segment`](@ref), [`KitePointGroup`](@ref), [`Pulley`](@ref) -""" -struct PointMassSystem - points::Vector{Point} - groups::Vector{KitePointGroup} - segments::Vector{Segment} - pulleys::Vector{Pulley} - tethers::Vector{Tether} - winches::Vector{Winch} - function PointMassSystem(points, groups, segments, pulleys, tethers, winches) - for (i, point) in enumerate(points) - @assert point.idx == i - end - for (i, group) in enumerate(groups) - @assert group.idx == i - end - for (i, segment) in enumerate(segments) - @assert segment.idx == i - end - for (i, pulley) in enumerate(pulleys) - @assert pulley.idx == i - end - for (i, tether) in enumerate(tethers) - @assert tether.idx == i - end - for (i, winch) in enumerate(winches) - @assert winch.idx == i - end - new(points, groups, segments, pulleys, tethers, winches) - end -end - """ mutable struct RamAirKite{S, V, P} <: AbstractKiteModel @@ -265,6 +53,10 @@ $(TYPEDFIELDS) e_z::V = zeros(S, 3) "Simplified system of the mtk model" sys::Union{ModelingToolkit.ODESystem, Nothing} = nothing + "Unsimplified system of the mtk model" + full_sys::Union{ModelingToolkit.ODESystem, Nothing} = nothing + "Linearization function of the mtk model" + lin_prob::Union{ModelingToolkit.LinearizationProblem, Nothing} = nothing "Velocity of the kite" vel_kite::V = zeros(S, 3) "Inertia around kite x y and z axis of the body frame" @@ -288,19 +80,26 @@ $(TYPEDFIELDS) iter::Int64 = 0 unknowns_vec::Vector{SimFloat} = zeros(SimFloat, 3) + defaults::Vector{Pair{Num, Real}} = Pair{Num, Real}[] + guesses::Vector{Pair{Num, Real}} = Pair{Num, Real}[] set_set_values::Function = () -> nothing set_measure::Function = () -> nothing set_vsm::Function = () -> nothing - set_unknowns::Function = () -> nothing - + set_unknowns::Function = () -> nothing + set_nonstiff::Function = () -> nothing + set_lin_vsm::Function = () -> nothing + set_lin_set_values::Function = () -> nothing + set_lin_unknowns::Function = () -> nothing + + get_vsm::Function = () -> nothing + get_set_values::Function = () -> nothing + get_unknowns::Function = () -> nothing get_state::Function = () -> nothing get_y::Function = () -> nothing prob::Union{OrdinaryDiffEqCore.ODEProblem, Nothing} = nothing - init_prob::Union{OrdinaryDiffEqCore.ODEProblem, Nothing} = nothing - integrator::Union{OrdinaryDiffEqCore.ODEIntegrator, Sundials.CVODEIntegrator, Nothing} = nothing - init_integrator::Union{OrdinaryDiffEqCore.ODEIntegrator, Sundials.CVODEIntegrator, Nothing} = nothing + integrator::Union{OrdinaryDiffEqCore.ODEIntegrator, Nothing} = nothing end function RamAirKite(set::Settings, aero::BodyAerodynamics, vsm_solver::VortexStepMethod.Solver, point_system::PointMassSystem) @@ -320,70 +119,90 @@ function RamAirKite(set::Settings, aero::BodyAerodynamics, vsm_solver::VortexSte return s end +function RamAirKite(set::Settings) + wing = RamAirWing(set; prn=false) + aero = BodyAerodynamics([wing]) + vsm_solver = Solver(aero; solver_type=NONLIN, atol=2e-8, rtol=2e-8) + point_system = PointMassSystem(set, wing) + return RamAirKite(set, aero, vsm_solver, point_system) +end + function update_sys_state!(ss::SysState, s::RamAirKite, zoom=1.0) - ss.time = s.t_0 - pos, acc, Q_b_w, elevation, azimuth, course, heading, e_x, tether_vel, twist, kite_vel = s.get_state(s.integrator) + isnothing(s.integrator) && throw(ArgumentError("run init_sim!(s) first")) + ss.time = s.integrator.t # Use integrator time + + # Get the extended state vector from the integrator + set_values, pos, acc_vec, Q_b_w, elevation, azimuth, course, heading_x, tether_length, tether_vel, winch_force, + twist_angle, kite_vel, aero_force_b, aero_moment_b, turn_rate, va_kite_b, wind_vec_gnd, wind_vel_kite = s.get_state(s.integrator) + P = length(s.point_system.points) for i in 1:P ss.X[i] = pos[1, i] * zoom ss.Y[i] = pos[2, i] * zoom ss.Z[i] = pos[3, i] * zoom end - # TODO - # ss.kite_acc .= kite_acc(s) - # ss.left_tether_vel = tether_vel[1] - # ss.right_tether_vel = tether_vel[2] - ss.acc = norm(acc) + + # --- Populate SysState fields --- + ss.acc = norm(acc_vec) # Use the norm of the kite's acceleration vector ss.orient .= Q_b_w + ss.turn_rates .= turn_rate ss.elevation = elevation ss.azimuth = azimuth - ss.force .= 0 - ss.heading = heading + + # Handle potential size mismatch for tether/winch related arrays + num_winches = length(s.point_system.winches) + ss.l_tether[1:num_winches] .= tether_length + ss.v_reelout[1:num_winches] .= tether_vel + ss.force[1:num_winches] .= winch_force + + # Depower and Steering from twist angles + num_groups = length(s.point_system.groups) + ss.twist_angles[1:num_groups] .= twist_angle + ss.depower = rad2deg(mean(twist_angle)) # Average twist for depower + ss.steering = rad2deg(twist_angle[end] - twist_angle[1]) + ss.heading = heading_x # Use heading_x from MTK model ss.course = course - ss.l_tether.= [s.tether_lengths; 0] - ss.v_reelout = [tether_vel; 0] - ss.depower = rad2deg(mean(twist)) - ss.steering = rad2deg(twist[end] - twist[1]) + # Apparent Wind and Aerodynamics + ss.v_app = norm(va_kite_b) + ss.v_wind_gnd .= wind_vec_gnd + ss.v_wind_kite .= wind_vel_kite + # Calculate AoA and Side Slip from apparent wind in body frame + # AoA: angle between v_app projected onto xz-plane and x-axis + # Side Slip: angle between v_app and the xz-plane + if ss.v_app > 1e-6 # Avoid division by zero + ss.AoA = atan(va_kite_b[3], va_kite_b[1]) + ss.side_slip = asin(va_kite_b[2] / ss.v_app) + else + ss.AoA = 0.0 + ss.side_slip = 0.0 + end + ss.aero_force_b .= aero_force_b + ss.aero_moment_b .= aero_moment_b ss.vel_kite .= kite_vel + # Calculate Roll, Pitch, Yaw from Quaternion + q = Q_b_w + # roll (x-axis rotation) + sinr_cosp = 2 * (q[1] * q[2] + q[3] * q[4]) + cosr_cosp = 1 - 2 * (q[2] * q[2] + q[3] * q[3]) + ss.roll = atan(sinr_cosp, cosr_cosp) + # pitch (y-axis rotation) + sinp = 2 * (q[1] * q[3] - q[4] * q[2]) + if abs(sinp) >= 1 + ss.pitch = copysign(pi / 2, sinp) # use 90 degrees if out of range + else + ss.pitch = asin(sinp) + end + # yaw (z-axis rotation) + siny_cosp = 2 * (q[1] * q[4] + q[2] * q[3]) + cosy_cosp = 1 - 2 * (q[3] * q[3] + q[4] * q[4]) + ss.yaw = atan(siny_cosp, cosy_cosp) + ss.set_torque[1:3] .= set_values nothing end -function SysState(s::RamAirKite, zoom=1.0) # TODO: add left and right lines, stop using getters and setters - isnothing(s.integrator) && throw(ArgumentError("run init_sim!(s) first")) - pos, acc, Q_b_w, elevation, azimuth, course, heading, e_x, tether_vel, twist, kite_vel = s.get_state(s.integrator) - P = length(s.point_system.points) - X = zeros(MVector{P, MyFloat}) - Y = zeros(MVector{P, MyFloat}) - Z = zeros(MVector{P, MyFloat}) - for i in 1:P - X[i] = pos[1, i] * zoom - Y[i] = pos[2, i] * zoom - Z[i] = pos[3, i] * zoom - end - - orient = MVector{4, Float32}(Q_b_w) # TODO: add Q_b_w - # forces = s.get_tether_force() # TODO: add tether force - forces = zeros(3) - t_sim = 0 - depower = rad2deg(mean(twist)) - steering = rad2deg(twist[end] - twist[1]) - ss = SysState{P}() - ss.time = s.t_0 - ss.t_sim = t_sim - ss.orient .= orient - ss.elevation = elevation - ss.azimuth = azimuth - ss.l_tether .= [s.tether_lengths; 0] - ss.v_reelout .= [tether_vel; 0] - ss.force .= [forces; 0] - ss.depower = depower - ss.steering = steering - ss.heading = heading - ss.course = course - ss.vel_kite .= kite_vel - ss.X = X - ss.Y = Y - ss.Z = Z +function SysState(s::RamAirKite, zoom=1.0) + ss = SysState{length(s.point_system.points)}() + update_sys_state!(ss, s, zoom) ss end @@ -415,35 +234,42 @@ and only update the state variables. Otherwise, it will create a new model from # Returns `Nothing` """ -function init_sim!(s::RamAirKite, measure::Measurement; prn=true, precompile=false) +function init_sim!(s::RamAirKite, measure::Measurement; + solver=ifelse(s.set.quasi_static, FBDF(nlsolve=OrdinaryDiffEqNonlinearSolve.NLNewton(relax=0.4, max_iter=1000)), FBDF()), + adaptive=true, prn=true, precompile=false, remake=false, reload=false, lin_outputs=Num[] +) function init(s, measure) - init_Q_b_w, R_b_w = measure_to_q(measure) - init_kite_pos = init!(s.point_system, s.set, R_b_w) + init_Q_b_w, R_b_w = measure_to_q(measure, s.wing.R_cad_body) + init!(s.point_system, s.set, R_b_w, init_Q_b_w) init_va_b = R_b_w' * [s.set.v_wind, 0., 0.] - sys, defaults, guesses = create_sys!(s, s.point_system, measure; init_Q_b_w, init_kite_pos, init_va_b) + inputs = create_sys!(s, s.point_system, measure; init_va_b) prn && @info "Simplifying the system" - prn && @time sys = structural_simplify(sys; additional_passes=[ModelingToolkit.IfLifting]) - !prn && (sys = structural_simplify(sys; additional_passes=[ModelingToolkit.IfLifting])) + prn ? (@time (sys, _) = structural_simplify(s.full_sys, (inputs, []))) : + ((sys, _) = structural_simplify(sys, (inputs, []))) s.sys = sys dt = SimFloat(1/s.set.sample_freq) if prn @info "Creating ODEProblem" - @time s.prob = ODEProblem(s.sys, defaults, (0.0, dt); guesses) + @time s.prob = ODEProblem(s.sys, s.defaults, (0.0, dt); s.guesses) else - s.prob = ODEProblem(s.sys, defaults, (0.0, dt); guesses) + s.prob = ODEProblem(s.sys, s.defaults, (0.0, dt); s.guesses) + end + if length(lin_outputs) > 0 + lin_fun, _ = linearization_function(s.full_sys, [inputs...], lin_outputs; op=s.defaults, guesses=s.guesses) + s.lin_prob = LinearizationProblem(lin_fun, 0.0) end - serialize(prob_path, s.prob) + serialize(prob_path, (s.prob, s.full_sys, s.lin_prob, s.defaults, s.guesses)) s.integrator = nothing + return nothing end prob_path = joinpath(KiteUtils.get_data_path(), get_prob_name(s.set; precompile)) - if !ispath(prob_path) + if !ispath(prob_path) || remake init(s, measure) end - try - reinit!(s, measure; precompile, prn) - catch e + success = reinit!(s, measure; solver, adaptive, precompile, reload, lin_outputs) + if !success rm(prob_path) @info "Rebuilding the system. This can take some minutes..." init(s, measure) @@ -452,10 +278,19 @@ function init_sim!(s::RamAirKite, measure::Measurement; prn=true, precompile=fal return nothing end + function init_sim!(::RamAirKite; prn=true) throw(ArgumentError("Use the function init_sim!(s::RamAirKite, measure::Measurement) instead.")) end +function linearize(s::RamAirKite; set_values=s.get_set_values(s.integrator)) + isnothing(s.lin_prob) && throw(ArgumentError("Run init_sim! with remake=true and lin_outputs=...")) + s.set_lin_vsm(s.lin_prob, s.get_vsm(s.integrator)) + s.set_lin_set_values(s.lin_prob, set_values) + s.set_lin_unknowns(s.lin_prob, s.get_unknowns(s.integrator)) + return solve(s.lin_prob) +end + """ reinit!(s::RamAirKite; prn=true, precompile=false) -> Nothing @@ -485,32 +320,37 @@ and only updates the state variables to match the current `measure`. # Throws - `ArgumentError`: If no serialized problem exists (run `init_sim!` first) """ -function reinit!(s::RamAirKite, measure::Measurement; prn=true, reload=true, precompile=false) - isnothing(s.point_system) && (s.point_system = PointMassSystem(s, s.wing)) - - init_Q_b_w, R_b_w = measure_to_q(measure) - init_kite_pos = init!(s.point_system, s.set, R_b_w) +function reinit!( + s::RamAirKite, measure::Measurement; + solver=ifelse(s.set.quasi_static, FBDF(nlsolve=OrdinaryDiffEqNonlinearSolve.NLNewton(relax=0.4, max_iter=1000)), FBDF()), + adaptive=true, + prn=true, + reload=true, + precompile=false, + lin_outputs=Num[] +) + isnothing(s.point_system) && throw(ArgumentError("PointMassSystem not defined")) + + init_Q_b_w, R_b_w = measure_to_q(measure, s.wing.R_cad_body) + init!(s.point_system, s.set, R_b_w, init_Q_b_w) - if isnothing(s.prob) + if isnothing(s.prob) || reload prob_path = joinpath(KiteUtils.get_data_path(), get_prob_name(s.set; precompile)) !ispath(prob_path) && throw(ArgumentError("$prob_path not found. Run init_sim!(s::RamAirKite) first.")) try - s.prob = deserialize(prob_path) + (s.prob, s.full_sys, s.lin_prob, s.defaults, s.guesses) = deserialize(prob_path) + length(lin_outputs) > 0 && isnothing(s.lin_prob) && throw(ArgumentError("lin_prob is nothing.")) catch e @warn "Failure to deserialize $prob_path !" - throw(e) + return false end end if isnothing(s.integrator) || !successful_retcode(s.integrator.sol) || reload t = @elapsed begin dt = SimFloat(1/s.set.sample_freq) - if s.set.quasi_static - solver = FBDF(nlsolve=OrdinaryDiffEqNonlinearSolve.NLNewton(relax=0.4, max_iter=1000)) - else - solver = FBDF() - end s.sys = s.prob.f.sys - s.integrator = OrdinaryDiffEqCore.init(s.prob, solver; dt, abstol=s.set.abs_tol, reltol=s.set.rel_tol, save_on=false, save_everystep=false) + s.integrator = OrdinaryDiffEqCore.init(s.prob, solver; + adaptive, dt, abstol=s.set.abs_tol, reltol=s.set.rel_tol, save_on=false, save_everystep=false) sym_vec = get_unknowns(s) s.unknowns_vec = zeros(SimFloat, length(sym_vec)) generate_getters!(s, sym_vec) @@ -518,30 +358,52 @@ function reinit!(s::RamAirKite, measure::Measurement; prn=true, reload=true, pre prn && @info "Initialized integrator in $t seconds" end - init_unknowns_vec!(s, s.point_system, s.unknowns_vec, init_Q_b_w, init_kite_pos) + init_unknowns_vec!(s, s.point_system, s.unknowns_vec) s.set_unknowns(s.integrator, s.unknowns_vec) - set_t!(s.integrator, 0.0) OrdinaryDiffEqCore.reinit!(s.integrator, s.integrator.u; reinit_dae=true) linearize_vsm!(s) - return nothing + return true end function generate_getters!(s, sym_vec) sys = s.sys c = collect - - set_set_values = setp(sys, sys.set_values) - set_measure = setp(sys, sys.measured_wind_dir_gnd) - set_vsm = setp(sys, c.([ + vsm_sym = c.([ sys.last_x, sys.last_y, sys.vsm_jac, - ])) - set_unknowns = setu(sys, sym_vec) + ]) - get_state = getu(sys, - [c(sys.pos), c(sys.acc), c(sys.Q_b_w), sys.elevation, sys.azimuth, sys.course, sys.heading_y, - c(sys.e_x), c(sys.tether_vel), c(sys.twist_angle), c(sys.kite_vel)] + set_set_values = setp(sys, sys.set_values) + set_measure = setp(sys, sys.measured_wind_dir_gnd) + set_vsm = setp(sys, vsm_sym) + set_unknowns = setu(sys, sym_vec) + set_nonstiff = setu(sys, get_nonstiff_unknowns(s)) + + get_vsm = getp(sys, vsm_sym) + get_set_values = getp(sys, sys.set_values) + get_unknowns = getu(sys, sym_vec) + get_state = getu(sys, + [c(sys.set_values), + c(sys.pos), # Particle positions + c(sys.acc), # Kite center acceleration vector (world frame) + c(sys.Q_b_w), # Orientation quaternion + sys.elevation, # Elevation angle + sys.azimuth, # Azimuth angle + sys.course, # Course angle + sys.heading_x, # Heading angle (based on body x-axis projection) + c(sys.tether_length), # Unstretched length per winch + c(sys.tether_vel), # Reeling velocity per winch + c(sys.winch_force), # Force at winch connection point per winch + c(sys.twist_angle), # Twist angle per group + c(sys.kite_vel), # Kite center velocity vector (world frame) + c(sys.aero_force_b), # Aerodynamic force (body frame) + c(sys.aero_moment_b), # Aerodynamic moment (body frame) + c(sys.turn_rate), # Angular velocity (body frame) + c(sys.va_kite_b), # Apparent wind velocity (body frame) + c(sys.wind_vec_gnd), # Ground wind vector (world frame) + c(sys.wind_vel_kite) # Wind vector at kite height (world frame) + ] ) get_y = getu(sys, sys.y) @@ -549,23 +411,37 @@ function generate_getters!(s, sym_vec) s.set_measure = (integ, val) -> set_measure(integ, val) s.set_vsm = (integ, val) -> set_vsm(integ, val) s.set_unknowns = (integ, val) -> set_unknowns(integ, val) - + s.set_nonstiff = (integ, val) -> set_nonstiff(integ, val) + + s.get_vsm = (integ) -> get_vsm(integ) + s.get_set_values = (integ) -> get_set_values(integ) + s.get_unknowns = (integ) -> get_unknowns(integ) s.get_state = (integ) -> get_state(integ) s.get_y = (integ) -> get_y(integ) + + if !isnothing(s.lin_prob) + set_lin_set_values = setp(s.lin_prob, sys.set_values) + set_lin_unknowns = setu(s.lin_prob, Initial.(sym_vec)) + set_lin_vsm = setp(s.lin_prob, vsm_sym) + + s.set_lin_set_values = (lin_prob, val) -> set_lin_set_values(lin_prob, val) + s.set_lin_unknowns = (lin_prob, val) -> set_lin_unknowns(lin_prob, val) + s.set_lin_vsm = (lin_prob, val) -> set_lin_vsm(lin_prob, val) + end nothing end -function linearize_vsm!(s::RamAirKite) - y = s.get_y(s.integrator) +function linearize_vsm!(s::RamAirKite, integ=s.integrator) + y = s.get_y(integ) jac, x = VortexStepMethod.linearize( s.vsm_solver, s.aero, y; - theta_idxs=1:4, - va_idxs=5:7, - omega_idxs=8:10, - moment_frac=s.bridle_fracs[s.point_system.groups[1].fixed_index]) - s.set_vsm(s.integrator, [x, y, jac]) + va_idxs=1:3, + omega_idxs=4:6, + theta_idxs=7:6+length(s.point_system.groups), + moment_frac=s.point_system.groups[1].moment_frac) + s.set_vsm(integ, [x, y, jac]) nothing end @@ -592,15 +468,12 @@ end function get_prob_name(set::Settings; precompile=false) suffix = "" - ver = "$(VERSION.major).$(VERSION.minor)_" + ver = "$(VERSION.major).$(VERSION.minor)" if precompile suffix = ".default" end - if set.quasi_static - return "prob_static_" * ver * string(set.segments) * "_seg.bin" * suffix - else - return "prob_dynamic_" * ver * string(set.segments) * "_seg.bin" * suffix - end + dynamics_type = ifelse(set.quasi_static, "static", "dynamic") + return "prob_$(ver)_$(set.physical_model)_$(dynamics_type)_$(set.segments)_seg.bin$suffix" end """ @@ -615,3 +488,111 @@ function calc_aoa(s::RamAirKite) return alpha_array[middle+1] end end + +function init_unknowns_vec!( + s::RamAirKite, + system::PointMassSystem, + vec::Vector{SimFloat} +) + !s.set.quasi_static && (length(vec) != length(s.integrator.u)) && + throw(ArgumentError("Unknowns of length $(length(s.integrator.u)) but vector provided of length $(length(vec))")) + + @unpack points, groups, segments, pulleys, winches, kite = system + vec_idx = 1 + + for point in points + if point.type == DYNAMIC + for i in 1:3 + vec[vec_idx] = point.pos_w[i] + vec_idx += 1 + end + for i in 1:3 + vec[vec_idx] = point.vel_w[i] + vec_idx += 1 + end + end + end + for pulley in pulleys + if pulley.type == DYNAMIC + vec[vec_idx] = pulley.length + vec_idx += 1 + vec[vec_idx] = pulley.vel + vec_idx += 1 + end + end + for group in groups + if group.type == DYNAMIC + vec[vec_idx] = group.twist + vec_idx += 1 + vec[vec_idx] = group.twist_vel + vec_idx += 1 + end + end + for winch in winches + vec[vec_idx] = winch.tether_length + vec_idx += 1 + vec[vec_idx] = winch.tether_vel + vec_idx += 1 + end + for i in 1:4 + vec[vec_idx] = kite.orient[i] + vec_idx += 1 + end + for i in 1:3 + vec[vec_idx] = kite.angular_vel[i] + vec_idx += 1 + end + for i in 1:3 + vec[vec_idx] = kite.pos[i] + vec_idx += 1 + end + for i in 1:3 + vec[vec_idx] = kite.vel[i] + vec_idx += 1 + end + + (vec_idx-1 != length(vec)) && + throw(ArgumentError("Unknowns vec is of length $(length(vec)) but the last index is $(vec_idx-1)")) + nothing +end + +function get_unknowns(s::RamAirKite) + vec = Num[] + @unpack points, groups, segments, pulleys, winches, kite = s.point_system + sys = s.sys + for point in points + for i in 1:3 + point.type == DYNAMIC && push!(vec, sys.pos[i, point.idx]) + end + for i in 1:3 # TODO: add speed to init + point.type == DYNAMIC && push!(vec, sys.vel[i, point.idx]) + end + end + for pulley in pulleys + pulley.type == DYNAMIC && push!(vec, sys.pulley_l0[pulley.idx]) + pulley.type == DYNAMIC && push!(vec, sys.pulley_vel[pulley.idx]) + end + vec = get_nonstiff_unknowns(s, vec) + !s.set.quasi_static && (length(vec) != length(s.integrator.u)) && + throw(ArgumentError("Integrator unknowns of length $(length(s.integrator.u)) should equal vec of length $(length(vec))")) + return vec +end + +function get_nonstiff_unknowns(s::RamAirKite, vec=Num[]) + @unpack points, groups, segments, pulleys, winches, kite = s.point_system + sys = s.sys + + for group in groups + group.type == DYNAMIC && push!(vec, sys.free_twist_angle[group.idx]) + group.type == DYNAMIC && push!(vec, sys.twist_ω[group.idx]) + end + for winch in winches + push!(vec, sys.tether_length[winch.idx]) + push!(vec, sys.tether_vel[winch.idx]) + end + [push!(vec, sys.Q_b_w[i]) for i in 1:4] + [push!(vec, sys.ω_b[i]) for i in 1:3] + [push!(vec, sys.kite_pos[i]) for i in 1:3] + [push!(vec, sys.kite_vel[i]) for i in 1:3] + return vec +end diff --git a/test/create_xz_file.jl b/test/create_xz_file.jl new file mode 100644 index 000000000..4be7ad911 --- /dev/null +++ b/test/create_xz_file.jl @@ -0,0 +1,49 @@ +# Copyright (c) 2024, 2025 Bart van de Lint, Uwe Fechner +# SPDX-License-Identifier: MIT + +using Timers +tic() +@info "Loading packages " + +using KiteModels, LinearAlgebra, Statistics + +toc() + +# Simulation parameters +dt = 0.05 +total_time = 10 # Longer simulation to see oscillations +vsm_interval = 3 +steps = Int(round(total_time / dt)) + +# Steering parameters +steering_freq = 1/2 # Hz - full left-right cycle frequency +steering_magnitude = 10.0 # Magnitude of steering input [Nm] + +# Initialize model +set = load_settings("system_ram.yaml") +set.segments = 3 +set_values = [-50, 0.0, 0.0] # Set values of the torques of the three winches. [Nm] +set.quasi_static = false +set.physical_model = "ram" + +@info "Creating wing, aero, vsm_solver, point_system and s:" +s = RamAirKite(set) +s.set.abs_tol = 1e-2 +s.set.rel_tol = 1e-2 +toc() + +# init_Q_b_w, R_b_w = KiteModels.measure_to_q(measure) +# init_kite_pos = init!(s.point_system, s.set, R_b_w, init_Q_b_w) +# plot(s.point_system, 0.0; zoom=false, front=true) + +measure = Measurement() + +# Initialize at elevation +s.point_system.winches[2].tether_length += 0.2 +s.point_system.winches[3].tether_length += 0.2 +measure.sphere_pos .= deg2rad.([70.0 70.0; 1.0 -1.0]) +KiteModels.init_sim!(s, measure; remake=false, reload=true) +sys = s.sys + +@info "System initialized at:" +toc() diff --git a/test/test_for_precompile.jl b/test/test_for_precompile.jl index 61eff460a..ebdfe3671 100644 --- a/test/test_for_precompile.jl +++ b/test/test_for_precompile.jl @@ -103,18 +103,14 @@ let simulate(integrator, 100, true) end if ! haskey(ENV, "NO_MTK") - using KiteModels, OrdinaryDiffEqCore, OrdinaryDiffEqBDF, OrdinaryDiffEqSDIRK, LinearAlgebra + using KiteModels, OrdinaryDiffEqCore, OrdinaryDiffEqBDF, LinearAlgebra using Base: summarysize update_settings() set = se("system_ram.yaml") set.segments = 2 set_values = [-50, -1.1, -1.1] - wing = RamAirWing(set) - aero = BodyAerodynamics([wing]) - vsm_solver = Solver(aero; solver_type=NONLIN, atol=1e-8, rtol=1e-8) - point_system = PointMassSystem(set, wing) - mtk_kite = RamAirKite(set, aero, vsm_solver, point_system) + mtk_kite = RamAirKite(set) measure = Measurement() KiteModels.init_sim!(mtk_kite, measure) logger = Logger(length(mtk_kite.point_system.points), 5) diff --git a/test/test_ram_air_kite.jl b/test/test_ram_air_kite.jl index 34e0417df..622c1440c 100644 --- a/test/test_ram_air_kite.jl +++ b/test/test_ram_air_kite.jl @@ -9,10 +9,6 @@ old_path = get_data_path() package_data_path = joinpath(dirname(dirname(pathof(KiteModels))), "data") temp_data_path = joinpath(tempdir(), "data") Base.Filesystem.cptree(package_data_path, temp_data_path; force=true) -for bin_file in ["ram_air_kite_body_info.bin", "ram_air_kite_foil_polar.bin"] - bin_path = joinpath(temp_data_path, bin_file) - isfile(bin_path) && rm(bin_path) -end set_data_path(temp_data_path) # Testing tolerance @@ -20,25 +16,24 @@ const TOL = 1e-5 const BUILD_SYS = true @testset verbose = true "RamAirKite MTK Model Tests" begin + # Initialize model set = se("system_ram.yaml") - wing = RamAirWing(joinpath(get_data_path(), "ram_air_kite_body.obj"), joinpath(get_data_path(), "ram_air_kite_foil.dat"); - mass=set.mass, crease_frac=0.82, align_to_principal=true) - aero = BodyAerodynamics([wing]) - vsm_solver = Solver(aero; solver_type=NONLIN, atol=1e-8, rtol=1e-8) - point_system = PointMassSystem(set, wing) + set.segments = 3 + set_values = [-50, 0.0, 0.0] # Set values of the torques of the three winches. [Nm] + set.quasi_static = false + set.physical_model = "ram" + + @info "Creating s:" + @time s = RamAirKite(set) + measure = Measurement() + s.set.abs_tol = 1e-2 + s.set.rel_tol = 1e-2 + + # Initialize at elevation + measure.sphere_pos .= deg2rad.([80.0 80.0; 1.0 -1.0]) - # Utility functions for setup - function create_test_model() - set.segments = 2 - set.abs_tol = TOL - set.rel_tol = TOL - VortexStepMethod.init!(aero; init_aero=false) - return RamAirKite(set, aero, vsm_solver, point_system) - end - @testset "Model Initialization Chain" begin - s = create_test_model() if BUILD_SYS # Delete existing problem file to force init! @info "Data path: $(get_data_path())" @@ -47,14 +42,14 @@ const BUILD_SYS = true @info "Removing existing serialized problem from $prob_path to test full initialization" rm(prob_path) end - + # 1. First time initialization - should create new model @info "Testing initial init! (should create new model)..." @time KiteModels.init_sim!(s, measure; prn=true) - + # Check that serialization worked @test isfile(prob_path) - + # Check initialization results @test !isnothing(s.integrator) @test !isnothing(s.sys) @@ -62,134 +57,141 @@ const BUILD_SYS = true end s.integrator = nothing s.sys = nothing - + # Keep references to first integrator and point system first_integrator_ptr = objectid(s.integrator) first_point_system_ptr = objectid(s.point_system) - + # 2. First reinit! - should load from serialized file @info "Testing first reinit! (should load serialized file)..." @time KiteModels.reinit!(s, measure; prn=true, reload=false) next_step!(s) - + # Check that it's a new integrator second_integrator_ptr = objectid(s.integrator) second_point_system_ptr = objectid(s.point_system) @test first_integrator_ptr != second_integrator_ptr @test first_point_system_ptr == second_point_system_ptr - + # 3. Second reinit! - should reuse existing integrator @info "Testing second reinit! (should reuse integrator)..." @time KiteModels.reinit!(s, measure; prn=true, reload=false) - + # This should create a new point system but reuse the existing integrator third_integrator_ptr = objectid(s.integrator) third_point_system_ptr = objectid(s.point_system) - @test second_integrator_ptr == third_integrator_ptr # Should be the same + @test second_integrator_ptr == third_integrator_ptr # Should be the same @test second_point_system_ptr == third_point_system_ptr - - # Get positions from various sources - pos_integrator, _, _, _, _, _, _, _, _, _, _ = s.get_state(s.integrator) + + # Get positions using SysState sys_state = KiteModels.SysState(s) - + # Check dimension consistency - @test size(pos_integrator, 2) == length(s.point_system.points) + # Note: pos_integrator is no longer directly fetched, comparing SysState to point_system @test length(sys_state.X) == length(s.point_system.points) - # Compare positions in different representations for (i, point) in enumerate(s.point_system.points) - # Points' world positions should match integrator positions + # Points' world positions should match SysState positions point_pos = point.pos_w - integ_pos = pos_integrator[:, i] sys_state_pos = [sys_state.X[i], sys_state.Y[i], sys_state.Z[i]] - - @test isapprox(norm(point_pos), norm(integ_pos), rtol=1e-2) - @test isapprox(norm(sys_state_pos), norm(integ_pos), rtol=1e-2) - + + # Use norm for comparison as exact vector match might be too strict due to float precision + @test isapprox(norm(point_pos), norm(sys_state_pos), rtol=1e-2) + # Positions should not be zero (except ground points) if point.type != KiteModels.WINCH # Skip ground points which might be at origin @test norm(point_pos) > 0.1 - @test norm(integ_pos) > 0.1 @test norm(sys_state_pos) > 0.1 end end end - + @testset "State Consistency" begin - s = create_test_model() KiteModels.reinit!(s, measure, prn=true, reload=false) - + sys_state_before = KiteModels.SysState(s) + # Check quaternion normalization - _, _, Q_b_w, elevation, azimuth, _, _, _, _, _, _ = s.get_state(s.integrator) - @test isapprox(norm(Q_b_w), 1.0, atol=TOL) - + @test isapprox(norm(s.integrator[s.sys.Q_b_w]), 1.0, atol=TOL) + # Check elevation matches measurement - @test isapprox(elevation, measure.elevation, atol=1e-2) - + @test isapprox(sys_state_before.elevation, measure.elevation, atol=1e-2) + # Change measurement and reinitialize old_elevation = measure.elevation measure.sphere_pos[1,:] .= deg2rad(85.0) KiteModels.reinit!(s, measure, prn=true, reload=false) - - # Get new state - _, _, _, elevation_new, _, _, _, _, _, _, _ = s.get_state(s.integrator) - + + # Get new state using SysState + sys_state_after = KiteModels.SysState(s) + # Verify state changed according to measurement - @test !isapprox(elevation_new, old_elevation, atol=1e-2) - @test isapprox(elevation_new, deg2rad(85.0), atol=1e-2) + @test !isapprox(sys_state_after.elevation, old_elevation, atol=1e-2) + @test isapprox(sys_state_after.elevation, deg2rad(85.0), atol=1e-2) end function test_step(s, d_set_values=zeros(3); dt=0.05, steps=5) + s.integrator.ps[s.sys.stabilize] = true + KiteModels.next_step!(s; dt=10.0) + s.integrator.ps[s.sys.stabilize] = false + @info "Stepping" for _ in 1:steps set_values = -s.set.drum_radius * s.integrator[s.sys.winch_force] + d_set_values KiteModels.next_step!(s, set_values; dt) + # Use SysState to get heading_x if needed, or directly from integrator if simpler + # sys_state_step = KiteModels.SysState(s) + # @show sys_state_step.heading_x # Example if heading_x is in SysState + @show s.integrator[s.sys.heading_x] # Keep direct access if simpler for this specific value end end - s = create_test_model() @testset "Simulation Step with SysState" begin # Basic step and time advancement test KiteModels.reinit!(s, measure; prn=true, reload=false) sys_state_before = KiteModels.SysState(s) - + # Run a simulation step with zero set values set_values = [0.0, 0.0, 0.0] dt = 1/s.set.sample_freq t, _ = KiteModels.next_step!(s, set_values; dt=dt) + # Update sys_state_before *after* the step to compare with the state *before* the loop KiteModels.update_sys_state!(sys_state_before, s) @test isapprox(t, dt, atol=TOL) - + # Run multiple steps num_steps = 10 total_time = 0.0 for _ in 1:num_steps step_time, _ = KiteModels.next_step!(s, set_values; dt=dt) - total_time += step_time + total_time += step_time # Accumulate time from next_step! return value end - sys_state_after = KiteModels.SysState(s) - @test any(abs.(sys_state_after.X .- sys_state_before.X) .> 0.1) - - @testset "Course Direction at Different Elevations" begin - # Test at 80 degrees elevation - measure.sphere_pos[1,:] .= deg2rad(50.0) - @test measure.elevation ≈ deg2rad(50.0) atol=1e-6 + sys_state_after = KiteModels.SysState(s) # Get state after the loop + # Compare state after loop with state after first step (stored in sys_state_before) + @test any(abs.(sys_state_after.X .- sys_state_before.X) .> 0.01) # Reduced tolerance slightly + + @testset "Course Direction at 60 Degrees Elevation" begin # Corrected description + # Initialize at 60 degrees elevation + measure.sphere_pos[1,:] .= deg2rad(60.0) + @test measure.elevation ≈ deg2rad(60.0) atol=1e-6 @test measure.azimuth ≈ 0.0 atol=1e-6 + KiteModels.reinit!(s, measure; prn=true) - @test s.integrator[s.sys.elevation] ≈ measure.elevation - @test s.integrator[s.sys.azimuth] ≈ measure.azimuth - test_step(s) - sys_state_50 = KiteModels.SysState(s) - @info "Course at 50 deg init elevation:" sys_state_50.course - @test sys_state_50.course ≈ 0.0 atol=π/4 - - # Test at 90 degrees elevation - measure.sphere_pos[1,:] .= deg2rad(89.0) - KiteModels.reinit!(s, measure; prn=true, reload=false) + + # Verify initial conditions using SysState + sys_state_init = KiteModels.SysState(s) + @test sys_state_init.elevation ≈ measure.elevation atol=1e-2 # Use relaxed tolerance consistent with other tests + @test sys_state_init.azimuth ≈ measure.azimuth atol=1e-2 + + # Run simulation steps test_step(s) - sys_state_89 = KiteModels.SysState(s) - @info "Course at 89 deg init elevation:" sys_state_89.course - @test abs(sys_state_89.course) ≈ π atol=π/4 + + # Check course direction using SysState + sys_state = KiteModels.SysState(s) + @info "Course at 60 deg elevation:" sys_state.course + + # At 60 degrees elevation, course should be roughly forward + @show sys_state.course + @test sys_state.course ≈ 0.0 atol=deg2rad(45.0) end # Utility function to calculate the signed angle difference between two angles @@ -200,36 +202,34 @@ const BUILD_SYS = true while diff < -π; diff += 2π; end return diff end - + @testset "Steering Response Using SysState" begin # Initialize model at moderate elevation - measure.sphere_pos[1,:] .= deg2rad(60.0) + measure.sphere_pos .= deg2rad.([70.0 70.0; 1.0 -1.0]) KiteModels.reinit!(s, measure; prn=true, reload=false) test_step(s) sys_state_initial = KiteModels.SysState(s) - + # steering right KiteModels.reinit!(s, measure; prn=true, reload=false) - test_step(s, [0, 0, -5]; steps=20) + test_step(s, [0, 10, -10]; steps=20) sys_state_right = KiteModels.SysState(s) - + # steering left KiteModels.reinit!(s, measure; prn=true, reload=false) - test_step(s, [0, -5, 0]; steps=20) + test_step(s, [0, -10, 10]; steps=20) sys_state_left = KiteModels.SysState(s) - + # Check steering values @info "Steering:" sys_state_right.steering sys_state_left.steering - @test sys_state_right.steering > 2.0 - @test sys_state_left.steering < -2.0 - @test abs(sys_state_right.steering) ≈ abs(sys_state_left.steering) atol=1.0 - + @test sys_state_right.steering > 3.0 + @test sys_state_left.steering < -3.0 + # Check heading changes right_heading_diff = angle_diff(sys_state_right.heading, sys_state_initial.heading) - @test right_heading_diff > 0.6 + @test right_heading_diff ≈ 0.6 atol=0.2 left_heading_diff = angle_diff(sys_state_left.heading, sys_state_initial.heading) - @test left_heading_diff < -0.5 - @test abs(right_heading_diff) ≈ abs(left_heading_diff) atol=0.3 + @test left_heading_diff ≈ -0.6 atol=0.2 end end end