Website in English

進化計算コンペティション2019
今年は風力発電用風車の設計問題です．

2019年12月14日(土) 9:00-12:00 兵庫県南あわじ市
主催：進化計算学会




コンペティション結果発表

【単目的最適化部門】優勝者　串田淳一（広島市立大学）
【多目的最適化部門】優勝者　野口隼（立命館大），原田智宏（首都大）

おめでとうございます！
参加者の皆様どうもありがとうございました．
コンペティションの結果については下記をご参照ください．

イントロダクション
コンペティション参加者のアルゴリズム紹介
風車の概念設計最適化問題
結果発表
結果の分析


お知らせ
・コンペティションの結果をアップロードしました．(2020/07/27)
・後処理プログラムにおいて，制約条件違反量が０であるときに制約条件を満たすと判定されるというミスがありましたので後処理プログラムを更新しました．なお，制約条件は実数であるため，この変更は結果には影響がありません．(2019/12/3)
・多目的最適化部門のハイパーボリューム算出のための正規化方法および参照点を修正しました．また，あわせて後処理プログラムを更新しました．(2019/11/27)
・制約条件の記述に一部誤りがありましたので修正しました．
（誤）「すべての制約条件はg(x)>=0の形で記述され，g(x)が0未満のとき制約条件を満足しないことを意味します．」
（正）「すべての制約条件はg(x)>0の形で記述され，g(x)が0より大きいとき制約条件を満足することを意味します．」(2019/11/27)
・ハイパーボリューム計算用コードについて，パレート解の数が４つ以下の場合にHypervolumeが正しく計算できないバグが確認されたので修正を行いました． (2019/11/6)
・ハイパーボリューム計算用コードについて，パレート解の数がちょうど3個の場合にHypervolumeが計算できないバグが確認されたので修正を行いました． (2019/10/29)
・Rで書かれたハイパーボリューム計算用コードでは計算時間がかかりすぎるという指摘がありましたので，C++版コードをアップロードいたしました．(2019/10/15)
・Q&Aを追加しました．(2019/10/4)
・同一名称の設計変数がなくなるように，【表１】の変数名に番号を振りました(2019/10/4)
・同一名称の制約条件がなくなるように，【表４】の制約条件名を変更しました(2019/10/4)
・Q&Aを追加しました．(2019/10/3)
・windturbine_SOP.pyおよびwindturbine_MOP.pyに結果には影響を与えない軽微なバグがあったので更新しました(2019/10/3)
・サンプルデータ(sample.zip)に間違いがあったので更新しました．(2019/10/3)
・制約条件の記述（表４）を修正しました．(2019/9/3)
・進化計算シンポジウム2019コンペティションのページを開設致しました．(2019/9/2)


コンペティション開催の趣旨
今年も進化計算の実問題での応用研究を促進するため最適化コンペティションを実施します． 例年通り，単目的設計最適化部門および多目的設計最適化部門の２部門です．
今回は再生可能エネルギーとして世界的に注目を集めている風力発電用風車の設計最適化問題を取り上げます．
風力発電用風車は発電量の増加や発電コストの低減のために，大型化が進んでいます． たとえば，NEDOプロジェクトの一環として（株）日立製作所が製作した5MW大型風力発電設備はローター直径126m，タワー高さ約90mあります． （参考）【国内最大級5MWの大型風力発電設備が完成】
一方でブレードの大型化に伴い翼端速度が速くなり，騒音などの新しい設計課題も生まれています．
今回取り扱う問題は実数の設計変数を３２個もつ，設計変数が比較的多い問題です．
また，制約条件も２２個あり，制約条件が厳しい問題でもあります．
多目的設計最適化部門では目的関数は５つ定義されており，皆さんお待ちかねの多数目的最適化問題となります．
上位入賞者には進化計算学会からの表彰もありますので，ぜひふるってご応募ください．

（参考）2017年度の第１回コンペティションのウェブサイト
（参考）2018年度の第２回コンペティションのウェブサイト

会場とアクセス
会場とアクセスについては進化計算シンポジウム2019のページをご覧ください． 進化計算シンポジウム2019の会場にて，初日(12/14)の午前中に実施します．（地元の方以外はおそらく前泊が必要になります）

主なスケジュール
コンペティション内容の公表および申込み受付開始：2019年9月2日（月）
申し込み締切： 2019年11月18日（月)（締め切りの延長はありません）
データ提出締切： 2019年12月2日（月)（締め切りの延長はありません）
コンペティション開催日：2019年12月14日（土）

申し込み方法
参加費は無料です．下記をご記入の上，2019年11月18日（月)までにメールでお申し込みをお願いします．
件名：進化計算シンポジウム2019コンペティション申し込み
本文：申込者の氏名，所属，メールアドレスと単目的部門・多目的部門のいずれへの申し込みかをご記入ください．
申込者は複数でもかまいません． 申込者が学生の場合は指導教員の名前も含めてください．
申込先： ec2019-competition@flab.isas.jaxa.jp

データ提出方法
2019年12月2日（月) までに下記に示すデータをメールでお送りください．
件名：進化計算シンポジウム2019コンペティションデータ送付（単目的部門 or 多目的部門）
申込先： ec2019-competition@flab.isas.jaxa.jp

暫定プログラム（変更の可能性があります）
9:00- 9:05 開催の趣旨説明　大山聖（宇宙航空研究開発機構）
9:05- 9:20 風車の概念設計最適化問題の説明 TBD
9:20-11:30 申込者からのプレゼン（１件あたり5～10分程度を予定）
11:30-12:00 結果の総括およびディスカッション　大山聖（宇宙航空研究開発機構）

表彰
各部門（単目的最適化および多目的最適化）の優勝者は進化計算学会から表彰があります．



ベンチマーク問題

今回のベンチマーク問題は風力発電用風車の設計最適化問題です．
【参考文献】Objectives and Constraints for Wind Turbine Optimization

設計変数および探査範囲を表１に示します．

【表１】設計変数および探査範囲

#下限設計変数名上限説明
11.0 [m]chord_sub[1]5.3 [m]ブレード翼弦長の翼幅方向分布[1]
21.0 [m]chord_sub[2]5.3 [m]ブレード翼弦長の翼幅方向分布[2]
31.0 [m]chord_sub[3]5.3 [m]ブレード翼弦長の翼幅方向分布[3]
41.0 [m]chord_sub[4]5.3 [m]ブレード翼弦長の翼幅方向分布[4]
50.1 [m]r_max_chord0.3 [m]ブレード最大翼弦長位置
6-5 [deg]theta_sub[1]30 [deg]ブレード取り付け角の翼幅方向分布[1]
7-5 [deg]theta_sub[2]30 [deg]ブレード取り付け角の翼幅方向分布[2]
8-5 [deg]theta_sub[3]30 [deg]ブレード取り付け角の翼幅方向分布[3]
9-5 [deg]theta_sub[4]30 [deg]ブレード取り付け角の翼幅方向分布[4]
100.005 [m]sparT[1]0.2 [m]ブレード翼桁厚さの翼幅方向分布[1]
110.005 [m]sparT[2]0.2 [m]ブレード翼桁厚さの翼幅方向分布[2]
120.005 [m]sparT[3]0.2 [m]ブレード翼桁厚さの翼幅方向分布[3]
130.005 [m]sparT[4]0.2 [m]ブレード翼桁厚さの翼幅方向分布[4]
140.005 [m]sparT[5]0.2 [m]ブレード翼桁厚さの翼幅方向分布[5]
150.005 [m]teT[1]0.2 [m]ブレード後縁補強材厚さの翼幅方向分布[1]
160.005 [m]teT[2]0.2 [m]ブレード後縁補強材厚さの翼幅方向分布[2]
170.005 [m]teT[3]0.2 [m]ブレード後縁補強材厚さの翼幅方向分布[3]
180.005 [m]teT[4]0.2 [m]ブレード後縁補強材厚さの翼幅方向分布[4]
190.005 [m]teT[5]0.2 [m]ブレード後縁補強材厚さの翼幅方向分布[5]
20-6.3 [m]precurve_sub[1]0.0 [m]ブレードの初期たわみの翼幅方向分布[1]
21-6.3 [m]precurve_sub[2]0.0 [m]ブレードの初期たわみの翼幅方向分布[2]
22-6.3 [m]precurve_sub[3]0.0 [m]ブレードの初期たわみの翼幅方向分布[3]
236.0tsr14.0 設計周速比(回転速度と風速の比)
246.0[rpm]maxOmega20.0[rpm]最大回転数
2550 [m]bladeLength80 [m]ブレード長さ
2620 [m]z_param70 [m]タワーくびれ位置
273.87 [m]tower_d[1]6.3 [m]タワー外径の高さ方向分布[1]
283.87 [m]tower_d[2]6.3 [m]タワー外径の高さ方向分布[2]
293.87 [m]tower_d[3]6.3 [m]タワー外径の高さ方向分布[3]
300.005 [m]t_param[1]0.1 [m]タワー厚さの高さ方向分布[1]
310.005 [m]t_param[2]0.1 [m]タワー厚さの高さ方向分布[2]
320.005 [m]t_param[3]0.1 [m]タワー厚さの高さ方向分布[3]


単目的設計最適化問題では発電コスト（＝（製造・建設コスト＋運用・保守コスト）/ 発電量）の最小化を設計目的とします．

【表２】設計目的（単目的最適化部門）
No. 目的関数 説明
1 発電コスト（＝（製造・建設コスト＋運用・保守コスト）/ 発電量）[$/kWh] 【最小化】発電コストを下げれば電気料金の引き下げにつながります．

多目的設計最適化問題は５つの設計目的を持ちます．

【表３】設計目的（多目的最適化部門）
No.目的関数 説明
1 年間発電量 [kWh/(year・turbine)] 【最大化】発電量を増やせば発電コスト削減や発電会社の利益の増加につながります．
2平均年間コスト（製造・建設コスト＋運用・保守コスト）[$/(year・turbine)] 【最小化】コストを下げれば発電コスト削減につながります．
3 タワー根元荷重 [Nm] 【最小化】タワー根元荷重を下げれば洋上に建設する場合の土台の建設コストが削減できます．
4 ブレード翼端速度 [m/s] 【最小化】ブレード翼端速度を下げれば空力騒音削減につながります．
5 疲労損傷度（ブレードおよびタワーの両方を考慮）[log scale] 【最小化】疲労損傷度を下げれば寿命が長くなり発電コスト削減につながります．


制約条件は表４に示します．

【表４】制約条件
No.制約条件 値 説明
1 ブレード最大たわみ量 [m] < ブレード・タワー間距離 / (安全率 1.485) [m] ブレードのたわみによるタワーとの衝突を防止
2 ブレードと地面との距離 [m] > 20 [m] 安全確保および地表付近の乱流を回避
3ブレードの固有振動数 [Hz] > ブレード通過周波数 * (安全率 1.1) [Hz] ブレードの共振を回避
4タワー固有振動数 [Hz] > ブレード回転周波数 * (安全率 1.1) [Hz] タワーの共振を回避
5 発電時のタワー最大応力 [Pa] < 4.5e+08 / (安全率 1.485) [Pa] 発電時の最大荷重での破損を防止
6 強風待機時のタワー最大応力 [Pa] < 4.5e+08 / (安全率 1.485) [Pa] 強風待機時の最大荷重での破損を防止
7 発電時のタワー全体座屈判定式 [-] < 1 (閾値) [-] 発電時の最大荷重での全体座屈による破損を防止
8 強風待機時のタワー全体座屈判定式 [-] < 1 (閾値) [-] 強風待機時の最大荷重での全体座屈による破損を防止
9 発電時のタワー局所座屈判定式 [-] < 1 (閾値) [-] 発電時の最大荷重での局所座屈による破損を防止
10 強風待機時のタワー局所座屈判定式 [-] < 1 (閾値) [-] 強風待機時の最大荷重での局所座屈による破損を防止
11 タワー疲労損傷度 [-] < 1 (閾値) [-] 20年以上のタワー寿命を確保
12タワー頂部直径 / タワー基部直径 [-] > 0.4 (閾値) [-] タワーの製造性を担保
13 タワー直径 / タワー厚さ [-] > 120 (閾値) [-] タワーの溶接性を担保
14 ブレード翼端速度 [m/s] < 80 [m/s] 風切り音を制限
15 ブレード翼桁ひずみ [-] > ブレード翼桁座屈ひずみ [-] ブレード翼桁の座屈による破損を防止
16 ブレード後縁ひずみ [-] > ブレード後縁座屈ひずみ [-] ブレード後縁補強材の座屈による破損を防止
17 ln(ブレード疲労損傷度) [-] < ln(1) = 0.0 (閾値) [-] 20年以上のブレード寿命を確保
18 ブレード翼桁圧縮ひずみ [-] > - (破断ひずみ 0.01) / (安全率 1.755) [-] ブレード翼桁の圧縮破壊を防止
19 ブレード翼桁引張ひずみ [-] < (破断ひずみ 0.01) / (安全率 1.755) [-] ブレード翼桁の引張破壊を防止
20 ブレード後縁圧縮ひずみ [-] > - (破断ひずみ 0.0025) / (安全率 1.755) [-] ブレード後縁補強材の圧縮破壊を防止
21 ブレード後縁引張ひずみ [-] < (破断ひずみ 0.0025) / (安全率 1.755) [-] ブレード後縁補強材の引張破壊を防止
22 発電量 [kWh] > 0.0 [kWh] 非物理的な解を省くための制約


コンペティションの条件

(1) 目的関数の評価にはNASAで開発されたOpenMDAOおよびNRELで開発されたWISDEMを用います．１つの設計候補の評価にかかる時間は約３秒です．
(2) 候補地点の評価回数の上限は一晩で最適化計算が終わることを目的として10,000回とします． １つの設計候補の評価にかかる時間は約３秒ですので，１最適化ケースの実施に約８時間かかることになります．
(3) 初期集団は設計空間全体に乱数で発生させる，実験計画法で計画的に配置する等，一般的な手法を用いてください．
(4) 評価基準
　・単目的最適化部門では発電コスト（=多目的最適化部門における目的関数2/目的関数1）が最も小さくなったアルゴリズムを最も優れているとします．
　　最終世代集団に含まれる解の最小値ではなく，最適化の過程で評価したすべての解の中で最も発電コストが小さい解の値とします．
　　発電コストが同じアルゴリズムが複数ある場合は，その発電コストを得るために要した設計評価回数が少ないアルゴリズムを最も優れているとします
　・多目的最適化部門ではハイパーボリューム(以下，HV)値が最も大きくなったアルゴリズムを最も優れているとします．
　　HV値は最終世代集団の非劣解のHV値ではなく，最適化の過程で評価したすべての解から作る集合の非劣解のHV値とします．
後述の評価モジュールは年間発電量ではなく，-1x年間発電量を出力しますので，すべての目的関数を最小化する問題になります．
(5) HV値の算出方法（2019/11/27更新）
ハイパーボリューム算出のための正規化はつぎのように行います．
$f_{normalized} = \frac{f - f^{min}}{f^{max} - f^{min}}$
$(f^{min}_1, f^{min}_2, f^{min}_3, f^{min}_4, f^{min}_5) = (-2.42e+07, 9.47e+05, 2.09e+07, 3.20e+01, -2.30e+00)$,
$(f^{max}_1, f^{max}_2, f^{max}_3, f^{max}_4, f^{max}_5) = (-1.96e+02, 2.35e+06, 1.89e+08, 8.00e+01, -2.50e-05)$.
HV値算出のための参照点は（1.5, 1.5, 1.5, 1.5, 1.5）とします． ただし更新された参照点を支配しない解については，ハイパーボリュームの算出には用いません． 詳細については下記のファイルを参照してください．
ハイパーボリュームの算出方法について

(6) 初期集団または最適化アルゴリズムの乱数を変えた21試行を実施し，その中央値となる発電コスト（単目的）またはHV値（多目的）で評価をします．


評価モジュールの実行方法

※単目的最適化/多目的最適化，共に同じ実行方法となります．

1. 下記のファイルのダウンロードと解凍をお願いします．
　[進化計算コンペティション2019インストールマニュアルver01.pdf] 初期インストールマニュアル
　[dummies.zip] ダミーモジュール
　[windturbine_SOP.py] 単目的最適化用評価モジュール
　[windturbine_MOP.py] 多目的最適化用評価モジュール
　[sample.zip] 計算結果例
2. 初期インストールマニュアルに従ってインストールをお願いします．
3. 評価モジュールを任意のフォルダに配置し，
python windturbine_sop.py (pop_vars_eval.txtのあるパス)
または
python windturbine_mop.py (pop_vars_eval.txtのあるパス)
で実行することができます．
4. 正常終了すると，pop_objs_eval.txtとpop_cons_eval.txtが生成されます．

pop_vars_eval.txt
評価モジュールが読み込む設計解データセット．1行1個体であり，タブで区切られた各列がそれぞれ設計変数を示します．
いずれも[0,1]で正規化しています．

pop_objs_eval.txt
評価モジュールが生成する目的関数のデータセット．1行1個体を示し，タブで区切られた各列がそれぞれ目的関数を示します．
行(個体)の順番はpop_vars_eval.txtと一致しています．
　各列は左から
　【単目的最適化】f1 = 発電コスト
　【多目的最適化】f1 = -年間発電量，f2 = 平均年間コスト，f3 = タワー根元荷重，f4 = ブレード翼端速度，f5 = 疲労損傷度
　です．多目的最適化におけるf1: 年間発電量は最大化すべき値のため，評価モジュール内部でマイナスをかけた値を出力しています．そのため全ての目的関数を最小化すればよい評価モジュールになっています．

pop_cons_eval.txt
評価モジュールが生成する制約条件のデータセット．1行1個体を示し，タブで区切られた各列がそれぞれ制約条件を示します．
行(個体)の順番はpop_vars_eval.txtと一致しています．
すべての制約条件はg(x)>0の形で記述され，g(x)が0より大きいとき制約条件を満足することを意味します．

実行方法で不明な点がありましたらec2019-competition@flab.isas.jaxa.jpまでお問い合わせください．


提出いただくデータ（単目的最適化部門）

(s1) 最適解の値 （21試行の中央値）（下記のフォーマットs1に従う）
(s2) 最適化の履歴（21試行の中央値）（下記のフォーマットs2に従う）
(s3) 最適解の値 （21試行の最良値）（下記のフォーマットs3に従う）
(s4) 最適化の履歴（21試行の最良値）（下記のフォーマットs4に従う）

【フォーマットs1】
- ファイルフォーマットはcsvファイル
- ファイル名は s_opt_***.csv 　（***は半角アルファベットで記載したグループ名）
- ベンチマーク問題の計算例にあるpareto.csvの書式に従って最適解の目的関数値，設計変数値，制約条件値などを記載してください．

【フォーマットs2】
- ファイルフォーマットはcsvファイル
- ファイル名は s_his_***.csv 　（***は半角アルファベットで記載したグループ名）
-　各列は左から，世代数，評価回数，その世代の最適解のf1=-発電コスト(，設計変数値，制約条件値),その世代で制約条件を満たした解の数，としてください．
　記載できない項目については-1を入れてください

【フォーマットs3】
- ファイルフォーマットはcsvファイル
- ファイル名は s_bst_***.csv 　（***は半角アルファベットで記載したグループ名）
- ベンチマーク問題の計算例にあるpareto.csvの書式に従って最適解の目的関数値，設計変数値，制約条件値などを記載してください．

【フォーマットs4】
- ファイルフォーマットはcsvファイル
- ファイル名は s_bst_his_***.csv 　（***は半角アルファベットで記載したグループ名）
-　各列は左から，世代数，評価回数，その世代の最適解のf1=-発電コスト(，設計変数値，制約条件値),その世代で制約条件を満たした解の数，としてください．
　記載できない項目については-1を入れてください



提出いただくデータ（多目的最適化部門）

(m1) HV値（21試行の中央値）
(m2) HV値（21試行の中央値）を算出する際に利用したすべての非劣解(劣解は除いてください）の設計変数，制約条件，目的関数など（フォーマットm2に従う）
(m3) 最適化（21試行の中央値）の履歴（フォーマットm3に従う）
(m4) HV値（21試行の最良値）
(m5) 最適化（21試行の最良値）の履歴（フォーマットm5に従う）

【フォーマットm2】
- ファイルフォーマットはcsvファイル
- ファイル名は m_prt_***.csv 　（***は半角アルファベットで記載したグループ名）
- ベンチマーク問題の計算例にあるpareto.csvの書式に従ってHV値の計算に利用した非劣解の目的関数値，設計変数値，制約条件値などを記載してください．
　記載できない項目については-1を入れてください

【フォーマットm3】
- ファイルフォーマットはcsvファイル
- ファイル名は m_his_***.csv 　（***は半角アルファベットで記載したグループ名）
- 各列は左から，世代数，評価回数，その世代集団のHV値，その世代までに評価したすべて解の集合のHV値，各世代ですべての制約条件を満たした実行可能解の数，としてください．
　記載できない項目については-1を入れてください

【フォーマットm5】
- ファイルフォーマットはcsvファイル
- ファイル名は m_bst_***.csv 　（***は半角アルファベットで記載したグループ名）
- 各列は左から，世代数，評価回数，その世代集団のHV値，その世代までに評価したすべて解の集合のHV値，各世代ですべての制約条件を満たした実行可能解の数，としてください．
　記載できない項目については-1を入れてください


後処理ツール

ベンチマーク問題が出力するファイルからの非劣解・HV値の算出のためのプログラム（C++コード）は下記からダウンロードしてください．
後処理プログラム（C++コード）(2019/12/03版）
※Rで書かれたハイパーボリューム計算用コードでは計算時間がかかりすぎるという指摘がありましたので，C++版コードに置換いたしました．

後処理にはご自分で開発されたプログラムを使っていただいてもかまいません．


その他

(1) コンペティション参加者には，コンペティション会場にて5～10分程度のプレゼンテーションをお願いします．
・利用したアルゴリズムの説明（工夫した点を含む）
・結果
について説明してください．
(2) 「結果の総括」において，みなさまから提出していただいたデータの比較を行いますのでご了承ください．
(3) 後日，解説記事等の執筆において，みなさまから提出していただいたデータを使わせていただく場合があります．もし問題がありましたらコンペティション前にお知らせください．


Q&A

Question 1: windturbine_SOP.pyやwindturbine_MOP.pyを実行すると下記の警告が出力されますが問題ないでしょうか？
UserWarning: rotor diameter may be modified in unexpected ways if tip precurve and precone are both nonzero
Answer: 問題ありませんので，そのままご使用ください．

Question 2: lcoe_se_assembly.pyを実行してもインストールマニュアルと同様の結果が出力されません．
Answer: pip freezeコマンドを実行し，OpenMDAOや各種Pythonパッケージのバージョンがインストールマニュアルで指定されたバージョンであるかご確認ください． また，pip freezeコマンドを実行し，WISDEMインストール時に導入したライブラリの欄に表示されるハッシュ値と，インストールマニュアルに記載された git checkoutコマンドで指定するハッシュ値が一致していることをご確認ください． バージョンやハッシュ値が一致していない場合，インストールマニュアルに従って再度インストール作業を行ってください．

Question 3:設計変数，目的関数，制約条件について有効桁数については考慮しないのですか？
Answer: 実問題については，設計変数，目的関数，制約条件の有効桁数について考慮すべきですが，今回のコンペティションでは有効桁数は考慮しないこととします．

Question 4: f5の正規化は-2.20から0.00の間で行われますが，f5の値が正にしかなりません．このまま計算をすすめてよいでしょうか？
Answer: f5は疲労損傷度のログスケールを取ったものです．f5が負の値の解は実行可能解，f5が正の値の解は実行不可能解です． 下に乱数で解を発生させた際のf5の値をプロットします． おおよそ10回に1回の実行可能解が現れますので，30回以上評価してすべて負になる場合はインストールがうまくいっていない可能性があります． その場合は，status_check.txtをec2019-competition@flab.isas.jaxa.jpにお送りください．





お問合せ先

ec2019-competition@flab.isas.jaxa.jp