View Javadoc
1   /*
2    * Licensed to the Hipparchus project under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The Hipparchus project licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      https://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  package org.hipparchus.optim.nonlinear.vector.constrained;
18  
19  import org.hipparchus.exception.LocalizedCoreFormats;
20  import org.hipparchus.exception.MathIllegalArgumentException;
21  import org.hipparchus.linear.RealMatrix;
22  import org.hipparchus.linear.RealVector;
23  import org.hipparchus.optim.LocalizedOptimFormats;
24  import org.hipparchus.optim.OptimizationData;
25  import org.hipparchus.optim.nonlinear.scalar.ObjectiveFunction;
26  import org.hipparchus.util.MathUtils;
27  
28  /**
29   * Abstract class for Sequential Quadratic Programming solvers
30   * @since 3.1
31   */
32  public abstract class AbstractSQPOptimizer extends ConstraintOptimizer {
33  
34      /** Algorithm settings. */
35      private SQPOption settings;
36  
37      /** Objective function. */
38      private TwiceDifferentiableFunction obj;
39  
40      /** Equality constraint (may be null). */
41      private EqualityConstraint eqConstraint;
42  
43      /** Inequality constraint (may be null). */
44      private InequalityConstraint iqConstraint;
45  
46      /** Simple constructor.
47       */
48      protected AbstractSQPOptimizer() {
49          this.settings = new SQPOption();
50      }
51  
52      /** Getter for settings.
53       * @return settings
54       */
55      public SQPOption getSettings() {
56          return settings;
57      }
58  
59      /** Getter for objective function.
60       * @return objective function
61       */
62      public TwiceDifferentiableFunction getObj() {
63          return obj;
64      }
65  
66      /** Getter for equality constraint.
67       * @return equality constraint
68       */
69      public EqualityConstraint getEqConstraint() {
70          return eqConstraint;
71      }
72  
73      /** Getter for inequality constraint.
74       * @return inequality constraint
75       */
76      public InequalityConstraint getIqConstraint() {
77          return iqConstraint;
78      }
79  
80      @Override
81      public LagrangeSolution optimize(OptimizationData... optData) {
82          return super.optimize(optData);
83      }
84  
85      @Override
86      protected void parseOptimizationData(OptimizationData... optData) {
87          super.parseOptimizationData(optData);
88          for (OptimizationData data : optData) {
89  
90              if (data instanceof ObjectiveFunction) {
91                  obj = (TwiceDifferentiableFunction) ((ObjectiveFunction) data).getObjectiveFunction();
92                  continue;
93              }
94  
95              if (data instanceof EqualityConstraint) {
96                  eqConstraint = (EqualityConstraint) data;
97                  continue;
98              }
99              if (data instanceof InequalityConstraint) {
100                 iqConstraint = (InequalityConstraint) data;
101                 continue;
102             }
103 
104             if (data instanceof SQPOption) {
105                 settings = (SQPOption) data;
106             }
107 
108         }
109 
110         // if we got here, convexObjective exists
111         int n = obj.dim();
112         if (eqConstraint != null) {
113             int nDual = eqConstraint.dimY();
114             if (nDual >= n) {
115                 throw new MathIllegalArgumentException(LocalizedOptimFormats.CONSTRAINTS_RANK, nDual, n);
116             }
117             int nTest = eqConstraint.dim();
118             if (nDual == 0) {
119                 throw new MathIllegalArgumentException(LocalizedCoreFormats.ZERO_NOT_ALLOWED);
120             }
121             MathUtils.checkDimension(nTest, n);
122         }
123 
124     }
125 
126     /**
127      * Compute Lagrangian gradient for variable X
128      *
129      * @param currentGrad current gradient
130      * @param jacobConstraint Jacobian
131      * @param x value of x
132      * @param y value of y
133      * @return Lagrangian
134      */
135     protected RealVector lagrangianGradX(final RealVector currentGrad, final RealMatrix jacobConstraint,
136                                          final RealVector x, final RealVector y) {
137 
138         int me = 0;
139         int mi;
140         RealVector partial = currentGrad.copy();
141         if (getEqConstraint() != null) {
142             me = getEqConstraint().dimY();
143 
144             RealVector ye = y.getSubVector(0, me);
145             RealMatrix jacobe = jacobConstraint.getSubMatrix(0, me - 1, 0, x.getDimension() - 1);
146 
147             RealVector firstTerm = jacobe.transpose().operate(ye);
148 
149             // partial = partial.subtract(firstTerm).add(jacobe.transpose().operate(ge).mapMultiply(rho));
150             partial = partial.subtract(firstTerm);
151         }
152 
153         if (getIqConstraint() != null) {
154             mi = getIqConstraint().dimY();
155 
156             RealVector yi = y.getSubVector(me, mi);
157             RealMatrix jacobi = jacobConstraint.getSubMatrix(me, me + mi - 1, 0, x.getDimension() - 1);
158 
159             RealVector firstTerm = jacobi.transpose().operate(yi);
160 
161             partial = partial.subtract(firstTerm);
162         }
163         return partial;
164     }
165 
166 }